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1.0.  INTRODUCTION 


Networks  of  queues  occur  frequently  In  diverse  applications.  In 
particular,  they  are  widely  used  In  studies  of  computer  and  communication 
system  performance  as  models  for  the  Interactions  among  system  resources. 
This  monograph  deals  with  mathematical  and  statistical  methods  for  discrete 
event  simulation  of  networks  of  queues.  The  emphasis  is  on  methods  for 
the  estimation  of  general  characteristics  of  "passage  times"  in  closed 
networks.  Informally,  a  passage  time  is  the  time  for  a  job  to  traverse 
a  portion  of  a  network.  Such  quantities,  calculated  as  random  sums  of 
queueing  times,  are  Important  in  computer  and  communication  system  models 
where  they  represent  job  response  times. 

By  simulation  we  mean  observation  of  the  behavior  of  a  stochastic 
system  of  Interest  by  artificial  sampling  on  a  digital  computer.  With 
discrete  event  simulation,  stochastic  changes  of  the  system  state  occur 
only  at  a  set  of  Increasing  time  points.  Simulation  is  a  tool  which  can 
be  used  to  study  complex  stochastic  systems  when  analytic  and/or  numerical 
techniques  do  not  suffice;  In  connection  with  the  study  of  complex  networks 
of  queues  encountered  in  applications,  this  is  often  the  case. 

When  simulating,  we  experiment  with  a  stochastic  system  and  observe 
Its  behavior.  During  the  simulation  we  measure  certain  quantities  In  the 
system  and,  using  statistical  techniques,  draw  Inferences  about 
characteristics  of  well  defined  random  variables.  The  most  obvious 
methodological  advantage  of  simulation  is  that  In  principle  it  is 
applicable  to  stochastic  systems  of  arbitrary  complexity.  In  practice, 


2 


however,  it  Is  often  a  decidedly  nontrivial  natter  to  obtain  from  a 
simulation  information  which  is  both  useful  and  accurate,  and  to  obtain 
it  in  an  efficient  manner.  These  difficulties  arise  primarily  from  the 
inherent  variability  in  a  stochastic  system,  and  it  is  necessary  to  seek 
theoretically  sound  and  computationally  efficient  methods  for  carrying 
out  the  simulation.  Apart  from  implementation  considerations.  Important 
concerns  for  simulation  relate  to  efficient  methoda  for  generation  of 
realizations  (sample  paths)  of  the  stochastic  system  under  study,  the 
design  of  simulation  experiments,  and  the  analysis  of  simulation  output. 

It  is  fundamental  for  simulation,  since  results  are  based  on  observation 
of  a  stochastic  system,  that  some  assessment  of  the  precision  of  results 
be  provided. 

Assessing  the  statistical  precision  of  a  point  estimate  requires 
careful  design  of  the  simulation  experiments  and  analysis  of  the  simulation 
output.  In  general,  the  desired  statistical  precision  takes  the  form  of 
a  confidence  Interval  for  the  quantity  of  interest.  Among  the  Issues  the 
simulator  must  face  are  the  initial  conditions  for  the  system  being 
simulated,  the  length  of  the  simulation  run,  the  number  of  replications 
of  the  experiments,  and  the  length  of  the  confidence  interval.  Over  the 
last  five  years,  there  has  been  Increased  attention  paid  to  these  issues, 
and  a  theory  of  simulation  analysis  (the  regenerative  method)  has  been 
developed  which,  when  applicable,  provides  some  measure  of  statistical 
precision.  The  regenerative  method,  which  is  based  on  limit  theorems 
developed  for  regenerative  stochastic  processes,  plays  a  key  role  in  our 
discussion  of  simulation  methods  for  passage  times  in  networks  of  queues. 
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Under  the  ueuel  queueing-theoretic  (independent  and  identically 
distributed  service  and  Interarrival  tine)  assumptions,  analyses  based  on 
a  "numbers-in-queue"  and  "stages-of-servlce"  state  vector  can  be  carried 
out.  Typically  it  is  necessary  to  assume  that  all  service  and  interarrival 
time  distributions  are  exponential  or  have  a  Cox-phase  (exponential  stage) 
representation.  Under  these  assumptions,  expressions  suitable  for 
numerical  evaluation  are  obtainable  for  queue  length  distributions.  Other 
measures  of  system  performance  (calculated  as  random  sums  of  queueing 
times)  Involve  the  times,  here  called  passage  times,  for  a  job  to  traverse 
a  portion  of  the  network.  Often  when  such  quantities  arise  in  computer 
and  communication  system  models,  they  represent  job  response  times.  In 
this  context,  characteristics  of  passage  times  other  than  expected  values 

(e.g.,  percentiles)  are  of  interest.  The  analyses  based  on  the 

* 

numbers-ln-queue,  stages-of-servlce  state  vector  yield  expected  values 
for  passage  times,  but  do  not  yield  other  passage  time  characteristics  of 
interest.  Moreover,  alternative  analyses  to  provide  these  measures  of 
the  variability  of  system  response  are  in  general  not  available,  and  it 
is  necessary  to  resort  to  simulation.  Although  the  usual  process  of 
numbers-in-queue  and  stages-of-service  is  a  regenerative  process  (in  fact 
a  Markov  chain)  under  the  probabilistic  assumptions  that  we  make  here, 
the  regenerative  method  cannot  be  applied  directly  to  this  process  to 
estimate  general  passage  time  characteristics.  This  is  essentially  because 
passage  times  are  not  totally  contained  within  cycles  of  the 
numbers-ln-queue,  stages-of-servlce  process. 
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The  organization  of  tha  presentation  is  as  follows.  This  initial 
section  provides  some  motivation  for  study  of  simulation  methods  for 
passage  times  in  networks  of  queues,  a  brief  overview  of  some  of  the 
methodological  considerations  for  simulation,  and  a  summary  of  the 
discussion  which  appears  in  subsequent  sections. 

The  estimation  methods  developed  here  for  passage  times  in  networks 
of  queues  use  the  regenerative  method  for  analysis  of  simulation  output. 
Based  on  a  single  simulation  run,  these  methods  provide  (strongly 
consistent)  point  estimates  and  (asymptotically)  valid  confidence  Intervals 
for  general  characteristics  of  limiting  passage  times.  Section  2  provides 
a  review  of  the  regenerative  method.  The  section  contains  a  brief 
discussion  of  the  underlying  theory  of  regenerative  stochastic  processes 
along  with  some  examples  of  regenerative  processes  in  networks  of  queues. 

Section  3  provides  a  specification  of  the  basic  class  of  closed 
networks  of  queues  with  which  we  deal,  and  the  probabilistic  assumptions 
therein.  Initially,  we  restrict  attention  to  networks  with  stochastically 
identical  jobs  and  give  a  state  vector  definition  based  on  a  linear  job 
stack.  The  section  also  contains  the  formal  definition  of  passage  time 
in  a  network  of  queues. 

The  notion  of  a  distinguished  "marked"  job  is  fundamental  to  the  method 
for  estimation  of  passage  time  characteristics  described  and  developed  in 
Section  4.  The  approach  is  to  consider  a  Markov  renewal  process  arising 
from  a  continuous  time  Markov  chain  defined  by  the  usual  number a -in-queue, 
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stages-of -service  state  vector  augmented  by  Information  sufficient  to 
track  the  marked  job.  We  arbitrarily  select  a  job  to  serve  as  the  marked 
job  and  measure  its  passage  times  during  the  simulation.  They  key  steps 
in  the  derivation  of  this  marked  job  method  are  identification  of  an 
appropriate  regenerative  process  in  discrete  time  and  development  of  a 
ratio  formula  from  which  point  estimates  and  confidence  intervals  can  be 
obtained  for  quantities  associated  with  the  limiting  passage  time. 

In  Section  5  we  consider  application  of  the  marked  job  method  to  two 
particular  closed  networks  of  queues ,  and  display  some  numerical  results. 
The  first  example  is  a  relatively  simple  network.  Despite  the  apparent 
structural  simplicity  of  this  network,  it  exhibits  the  essence  of  the 
passage  time  simulation  problem.  The  second  and  more  complex  network 
arises  as  a  model  for  a  computer  data  base  management  system.  This  model 
illustrates  the  representation  of  complex  congestion  phenomena  in  the 
framework  of  Section  3. 

The  extension  of  the  marked  job  method  to  certain  finite  capacity  open 
networks  of  queues  is  the  subject  of  Section  6.  Particular  stochastic 
point  processes  associated  with  a  Markov  renewal  process  generate  arrivals 
to  the  networks,  and  there  are  two  formulations  of  the  finite  capacity 
constraint.  The  network  structure  we  permit  is  essentially  the  same  as 
that  described  in  Section  3  except  that  here  the  networks  are  open.  To 
estimate  passage  times  in  these  networks,  we  track  an  appropriate  sequence 
of  typical  jobs,  based  on  the  idea  of  a  marked  job.  These  are  to  be 
typical  jobs  in  the  sense  that  the  sequence  of  passage  times  for  the  marked 
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jobs  should  converge  in  distribution  to  the  same  random  variable  as  do 
the  passage  times  for  all  the  jobs.  It  is  necessary  to  take  some  care  to 
ensure  chat  this  is  the  case. 

Slightly  restricting  the  definition  of  passage  time  given  in  Section  3, 
we  develop  in  the  next  section  a  new  and  somewhat  simpler  stochastic 
setting  than  that  of  Section  4  for  the  marked  job  method.  Here  the  starts 
of  passage  times  for  the  marked  job  are  the  successive  entrances  of  a 
particular  continuous  time  Markov  chain  to  a  fixed  set  of  states,  and  the 
terminations  of  such  passage  times  are  the  successive  entrances  to  another 
fixed  set  of  states.  This  f ovulation,  in  terms  of  hitting  times  to  fixed 
sets  of  states,  is  the  basis  for  the  passage  time  simulation  method 
developed  in  the  next  section. 

The  marked  job  method  of  Section  5  is  applicable  to  passage  times  in 
the  general  sense,  i.e.,  whether  or  not  the  passage  time  is  a  complete 
circuit.  For  those  passage  times  (termed  "response  times")  which  are 
complete  circuits  in  a  closed  network,  simulation  using  a  marked  job 
appears  to  be  the  only  method  available  for  obtaining  confidence  intervals 
from  a  single  simulation  run.  It  is  inherent  in  the  marked  job  method 
that  only  the  passage  times  observed  for  the  marked  job  enter  into  the 
construction  of  point  and  interval  estimates,  and  we  would  expect  some 
loss  of  efficiency  as  the  price  for  obtaining  confidence  intervals.  In 
Section  8,  we  concentrate  on  passage  times  through  a  subnetwork  of  a  given 
network  of  queues ,  and  develop  the  decomposition  method .  With  this 
estimation  method,  passage  times  observed  for  all  the  jobs  during  a  single 
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simulation  run  enter  Into  the  construction  o£  point  and  Interval  estimates. 
The  baels  for  this  method  Is  the  observation  that  the  successive  entrances 
to  a  fixed  set  of  states  of  an  appropriate  continuous  time  stochastic 
process  serve  to  decompose  the  sequence  of  passage  times  for  all  the  jobs 
into  Independent  and  identically  distributed  blocks. 

Section  9  deals  with  the  statistical  efficiency  of  the  marked  job  and 
decomposition  methods.  We  consider  the  calculation  of  theoretical  values 
for  variance  constants  entering  into  central  limit  theorems  used  to  obtain 
confidence  intervals  for  mean  passage  times.  The  results  of  this  section 
provide  a  firm  basis  for  comparing  the  efficiency  of  the  two  methods  where 
both  apply.  They  can  also  be  used  to  give  some  idea  of  the  efficiency  of 
the  marked  job  method  in  the  case  of  response  times,  where  it  is  our  only 
means  of  obtaining  confidence  intervals  from  a  single  simulation  run. 

The  estimation  of  passage  times  in  closed  networks  of  queues  with 
multiple  job  types  is  the  subject  of  Section  10.  Here  the  type  of  a  job 
may  influence  its  routing  through  the  network  as  well  as  its  service 
requirements  at  each  center.  Using  the  stochastic  setting  of  Section  7 
for  the  marked  job  method,  we  mark  one  job  of  each  type.  By  tracking 
these  marked  jobs  through  the  network,  we  are  able  to  obtain  point  and 
interval  estimates  for  a  variety  of  measures  of  the  variability  of  system 
response  over  the  several  job  types. 


The  final  section  considers  some  aspects  of  uniform  and  nonuniform 
random  number  generation  pertinent  to  implementation  of  a  passage  time 
simulation.  We  also  discuss  the  use  of  random  number  streams  and  the 
generation  of  state  vector  processes. 
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2.0.  SIMULATION  OF  REGENERATIVE  PROCESSES 

When  the  output  of  a  discrete  event  digital  simulation  is  a  stochastic 
process  g-{X(t):tiO)  that  approaches  a  "steady  state"  which  is  of  interest, 
it  may  be  possible  to  characterize  the  stochastic  structure  of  the  process 
and  to  use  this  structure  in  carrying  out  the  simulation.  In  such  cases, 
mathematical  results  on  the  stochastic  structure  of  the  process  g  form 
the  basis  for  both  the  design  of  simulation  experiments  and  the  analysis 
of  simulation  output.  For  particular  stochastic  processes  known  as 
regenerative  processes,  Fishman  (1973)  and  Crane  and  Iglehart  (1975a)  have 
provided  a  theory  of  simulation  analysis  called  the  regenerative  method. 
This  theory  has  been  developed  in  subsequent  papers  including  Crane  and 
Iglehart  (1975b),  Iglehart  (1975),  (1976),  Hordljk,  Iglehart  and 
Schassberger  (1976) ,  and  Lavenberg  and  Sauer  (1977) ;  see  Crane  and  Lemolne 
(1977)  and  Iglehart  (1978)  for  an  introduction  to  and  a  detailed  review 
of  the  regenerative  method.  Regenerative  simulation  underlies  the 
estimation  methods  described  in  subsequent  sections  for  passage  times  in 
networks  of  queues.  In  this  section  we  discuss  regenerative  processes 
and  review  the  regenerative  method. 

Heurlstically,  a  regenerative  stochastic  process  is  a  process  having 
the  property  that  there  exist  random  time  points  at  which  the  process 
probabilistically  restarts.  Typically,  these  time  points  at  which  the 
process  probabilistically  starts  afresh,  referred  to  as  regeneration  points 
or  regeneration  times,  are  returns  to  a  fixed  state  of  the  process.  The 
essential  idea  of  a  regenerative  process  is  that  between  any  two  successive 
regeneration  points,  the  evolution  of  the  process  is  a  probabilistic 
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replica,  of  the  process  between  any  ocher  such  pair  of  regeneration  points. 
For  many  stochastic  models  In  which  It  is  possible  to  Identify  a  sequence 
of  regeneration  points,  this  discovery  Is  a  key  to  an  analytic/numerical 
solution  which  yields  expressions  for  quantities  of  interest.  It  is  often 
the  case,  however,  that  even  though  a  sequence  of  regeneration  points 
exists,  it  is  nevertheless  not  possible  to  obtain  an  analytic/nuaerlcal 
solution  because  of  severe  computational  difficulties;  we  would  then 
consider  simulation,  and  it  is  this  situation  with  which  we  deal  here. 

The  regenerative  process  structure,  in  the  presence  of  certain 
regularity  conditions,  guarantees  the  existence  of  a  "steady  state"  for 
the  process,  i.e.,  that  these  exists  a  random  variable  X  such  that 

lim  P{X(t)Sx}  -  P{X£x}  , 

for  all  x  at  which  the  right  hand  side  of  this  equation  is  continuous. 

This  type  of  convergence  (in  distribution)  is  known  as  weak  convergence 
and  we  denote  it  by  "X(t)—>X  as  t-*»."  Furthermore,  the  regenerative 
structure  ensures  that  the  "steady  state"  X  of  the  process  is  determined 
(as  a  ratio  of  expected  values)  by  the  behavior  of  the  process  between 
two  successive  regeneration  points.  There  is  an  important  implication  of 
these  mathematical  results  for  the  simulation  of  regeneration  processes. 

A  strongly  consistent  point  estimate  and  asymptotically  valid  confidence 
interval  for  the  expected  value  of  a  general  (measurable)  function  of  the 
steady  state  X  can  be  obtained  by  observation  (in  cycles  of  random  length 
defined  by  the  regeneration  points)  of  a  finite  portion  of  a  single 
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realization  of  the  process  $.  When  simulating  (say,  for  a  fixed  number 
of  cycles) ,  we  measure  appropriate  quantities  defined  totally  within  the 
individual  cycles  and  compute  sample  means  over  the  cycles. 

Where  applicable,  the  regenerative  method  has  great  appeal  because  it 
provides  both  point  and  interval  estimates  having  desirable  properties. 
There  are,  however,  other  considerations.  The  classical  alternative  for 
estimation  of  "steady  state"  quantities  would  entail  selecting  an  initial 
state  for  the  process,  running  the  simulation  for  an  initial  period  of 
time  (and  discarding  this  "initial  transient") ,  and  then  observing  the 
process  ("in  steady  state")  for  an  additional  period  of  time  from  which 
point  estimates  are  obtained.  In  general,  no  confidence  Interval  is 
available,  nor  is  there  any  guidance  on  the  selection  of  the  initial  state. 
Moreover,  the  determination  of  appropriate  initial  and  additional  periods 
of  time  is  often  nontrivial  and  likely  to  require  sophisticated  statistical 
techniques.  There  are  similar  problems  with  simulation  methods  based  on 
multiple  replications.  With  the  regenerative  method,  these  difficulties 
to  a  large  extent  are  avoidable. 

2.1.  Definition  of  Regenerative  Process 

A  regenerative  process  in  continuous  time  can  be  defined  in  terms  of 
the  pasting  together  of  so-called  "tours";  see  Smith  (1958)  and  Miller 
(1972).  The  formal  definition  of  a  regenerative  process  that  we  give  is 
equivalent  to  these  and  also  to  the  definition  of  (Jinlar  (1975a),  p.  298. 

We  require  the  notion  of  a  renewal  process  and  that  of  a  stopping  time 
for  a  stochastic  process. 
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A  sequence  of  random  variables  X"^TQ:naO}  is  a  renewal  process  provided 

Chat  T.-O  and  T  -T  .  (nil)  are  independent,  Identically  distributed 
u  n  n— x 

(l.l.d.)  positive  random  variables.  We  always  assume  that  X  Is  persistent, 

i.e.,  that  P{T  -T  <»}■!.  a  random  variable  T  taking  values  in  [0,+») 
n  n-i 

is  a  stopping  time  for  a  stochastic  process  X  provided  that  for  every 
finite  tkO,  the  occurrence  or  nonoccurrence  of  the  event  {T£t}  can  be 
determined  from  the  history  {X(u):u£t}  of  the  process  up  to  time  t;  see 
Clnlar  (1975a),  p.  239  for  a  discussion  of  stopping  times. 

(2.1.1)  DEFINITION.  The  real  (possibly  vector-valued)  stochastic  process 
X*{ X ( t ) :tiO}  is  a  regenerative  process  if 

(1)  there  exists  a  sequence  of  stopping  times  ^-{B^rn^O}  such  that 
jg  is  a  renewal  process;  and 

(ii)  for  every  sequence  of  times  0<t^<t2<. . .<tm  (nfcl)  and  n^O,  the 

random  vectors  {X(t. ) , . . . ,X(t  )}  and  (X(B+t,) . X(B+t  )}  have 

i  ni  n  x  n  id 

the  same  distribution,  and  the  processes  {X(t):t<8n)  and 
{X(BQ+t) :t20)  are  Independent. 

The  points  of  £  are  the  regeneration  points  for  the  process  X  and  we  refer 

to  the  interval  [8  . ,8  )  as  the  nth  cycle  of  the  regenerative  process. 

n-JL  n 

The  definition  of  a  discrete  time  regenerative  process  is  similar;  see 
the  discussion  of  recurrent  events  in  Feller  (1968),  Ch.  XIII. 

It  is  straightforward  to  check  that  irreducible  and  positive  recurrent 
(continuous  or  discrete  time)  Markov  chains  and  semi-Markov  processes 
having  finite  (or  countable)  state  space  are  regenerative  processes.  It 
can  also  be  shown  that  in  a  single  server  queueing  system  in  which  the 
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sequence  of  successive  lnterarrlval  and  service  times  are  positive,  l.l.d. 

random  vectors  (having  finite  means),  the  processes  of  the  number  of  jobs 

Q(t)  in  the  system  at  time  t,  the  waiting  time  W  of  the  nth  job,  and  the 

n 

virtual  waiting  time  V(t)  at  time  t  are  all  regenerative  processes, 
provided  that  the  traffic  intensity  is  less  than  1.  For  the  process 
{WQ:n20},  the  regeneration  points  are  the  indices  of  jobs  having  zero 
waiting  time.  For  the  process  {Q(t):t*0}  as  well  as  the  process 
{V(t):t20},  the  regeneration  points  are  the  start  of  a  period  during  which 
the  server  is  busy. 

In  developing  the  regenerative  method,  it  is  necessary  to  distinguish 

two  cases.  For  the  renewal  process  £■{ 0Q :n20}  in  the  definition  of  the 

regenerative  process  J,  we  let  an“8n-6n_1  (n2l)  and  denote  by  F  the  distribution 

function  of  a  .  The  random  variable  a.  (or  distribution  function  F)  is 
n  x 

said  to  be  periodic  with  period  X>0  if,  with  probability  one,  assumes 

values  in  the  set  {0,X,2X . }  and  X  is  the  largest  such  number.  If 

there  is  no  such  X,  then  (or  F)  is  said  to  be  aperiodic . 

In  the  aperiodic  case,  in  order  for  a  regenerative  process  to  have  a 

limiting  distribution,  it  is  necessary  either  to  impose  regularity 

conditions  on  the  sample  paths  of  the  process,  or  to  place  restrictions 

on  the  distribution  function  of  the  time  between  regeneration  points.  To 

be  more  specific,  we  first  define  an  appropriate  class  of  distribution 

functions.  Let  F  be  the  n-fold  convolution  of  the  distribution  function  F, 
n 

and  define  d  to  be  the  set  of  all  distribution  functions  F  such  that  for 
some  n2l,  Fq  has  an  absolutely  continuous  component  (l.e.,  has  a  density 
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function  on  son*  Interval) .  It  le  convenient  to  write  when  the 

dlatrlbutlon  function  F  of  ie  en  element  of  Jtd .  We  would  expect  the 
aperiodic  distrlbutlona  F  arising  in  applications  to  be  elements  of  the 
set  .ri.  With  respect  to  appropriate  regularity  conditions  on  the  sample 
paths  of  the  process,  we  restrict  attention  to  processes  X«{X(t) :t20} 
having  right  continuous  sample  paths  and  limits  from  the  left,  l.e.,  for 
t20. 


X(t)  -  lim  X(u) 
u+t 

and  for  all  t>0, 

X(t-)  -  lim  X(u)  exists  , 
u+t 

with  probability  one.  For  such  a  k -dimensional  stochastic  process  we 
write  •  With  these  definitions,  we  can  state  the  basic  limit 

theorem  for  regenerative  processes. 

(2.1.2)  THEOREM.  Assume  that  is  aperiodic  with  E{a1}<».  If  either 

geDjJO,®)  or  then  X(t)->X  as  t-*». 

There  is  a  corresponding  result  for  the  periodic  case.  The  proof  of  this 

theorem  (Miller  (1972))  involves  an  application  of  the  key  renewal  theorem 

and  is  somewhat  involved  technically. 

Now  suppose  that  is  aperiodic  end  that  for  a  real-valued 
(measurable)  function  f  having  domain  E,  the  state  space  of  the  process 
J,  the  quantity  of  interest  is 


»•  >* 


15 


r(f)  -  E{f (X) }  . 


When  the  state  space  E  is  not  discrete,  we  must  also  assume  that  the  set 
D(f)  of  discontinuities  of  f  is  such  that  P{X«D(f)}*0;  if  E  is  discrete,  we 
can  choose  f  arbitrarily.  We  always  assume  that  with  probability  one  the 
process  f(X)  is  integrable  over  a  finite  interval,  and  for  nsl  define 


Yn(f) 


f (X(u))du  . 


In  the  case  of  a  discrete  time  regenerative  process,  for  n^l  we  define 


Vf>  ■  £  £<V  • 

k-e»-i 

Theorems  (2.1.3)  and  (2.1.4)  deal  with  the  structure  of  regenerative 
processes.  These  results  form  the  basis  for  the  regenerative  method. 


(2.1.3)  THEOREM.  The  sequence  { (YQ(f) ,aQ) :n^l)  consists  of  Independently 
and  Identically  distributed  random  vectors. 


This  follows  directly  from  the  definition  of  a  regenerative  process.  The  next 
result  (cf.  Crane  and  Iglehart  (1975a),  Appendix)  provides  a  ratio  formula 
for  the  quantity  r(f). 


(2.1.4)  THEOREM.  Assume  that  ot^  is  aperiodic  with  and  that 

E{  | f  (X)  |  }<».  If  either  fCDsD^o,-)  or  then 

E{f(X)}  -  E{Y1(f))/E{a1>  . 
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There  le  an  analogous  ratio  formula  when  ol  la  periodic.  Note  that  If 
the  state  space  E  is  discrete,  the  condition  f (g)eD^[Ot«)  always  holds. 

We  now  indicate  how  to  obtain  a  point  estimate  and  confidence  Interval 
for  the  quantity  r(f)»E{f(X)}  from  a  sample  path  of  the  process  g.  In 
this  discussion,  we  assume  that  the  regenerative  process  X  and  the  function 
f  are  such  that  the  ratio  formula  for  r(f)  holds.  For  kkl,  let 

Zk(f)  -  Yk(f)  -  r(f)ok  (2.1.5) 

and  denote  the  variance  of  Z^(f)  by  a^«var{Z^(f)} .  Note  that  {Zk(f) :kkl} 
consists  of  1.1. d.  random  variables,  that  Z^Cf)  la  completely  determined 
by  the  kth  cycle  of  the  regenerative  process  g,  and  that  E(Zk(f)}-0. 

Writing 

EUZ^f))2}  -  E{[(Y1(f)-E{Y1(f)»-r(f)(a1-E{a1})]2}  , 
it  follows  that 

a2  -  var{Y1(f)>  -  2r(f)  coviY^f)  ,0^}  +  (r(f))2  var{a1>  .  (2.1.6) 

2  2 

We  require  the  further  assumption  that  0<o  <°°.  The  case  a  "0  Is 
2 

degenerate,  and  a  «=°  for  most  finite  state  processes.  (In  some  queueing 

systems,  however,  additional  finite  higher  moment  conditions  on  service 

2 

and  lnterarrlval  times  are  needed  to  ensure  a  «*>.)  For  fixed  n,  we  use 

—  —  2 
the  quantities  Yn(f),  <*n,  a^t  ®22  and  ®12  t0  estio*te  a  *  obtained  from 

»®j_)  *  •  •  •  *  (^n(f)  »  these  are  the  usual  unbiased  estimators  of 

E{Y1(f)},  £{0^},  var{Y1(f)},  varfo^}  and  cov{Y1(f) ,0^},  respectively.  As 
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*  consequence  of  the  strong  law  of  larga  numbers  for  l.i.d.  sequences  of 
random  variables,  as  ir*®  the  point  estimates 


f  (f)  -  Y  (f )/a 
n  n  n 


•a  '  <*u-2V«-12+<V»>2*22)1/2 


converge  with  probability  one  to  r(f)  and  c,  respectively;  thus,  by 
definition  fQ(f)  and  sq  are  strongly  consistent  estimates. 


The  basis  for  the  construction  of  an  asymptotically  (ir*®)  valid 

confidence  interval  for  r(f)  is  a  particular  central  limit  theorem 
2 

(c.l.t.):  if  0<ct  <“,  then  as  tr*®, 


(2.1.7) 


n1/2{?  (fJ-rCfJJ/fa/ECa.)]  N(0,1)  , 
n  i 


where  N(0,1)  is  the  standardized  (mean  0,  variance  1)  normal  random 
variable. 


To  derive  this  result  and  similar  c.l.t.'s  later,  we  need  two  lemmas 

on  weak  convergence.  Let  (X  :n2l}  and  {Y  :n^l)  be  two  sequences  of  random 

n  n 

vectors  such  that  X  and  Y  are  defined  on  a  common  probability  space  for 
n  n 

k  i  k  r 

all  n  and  have  ranges  R  and  R  ,  where  R  [respectively  R  J  is 

k-d linens iona  1  [respectively  1-dimensional]  Euclidean  space.  Let  c  denote 

k  1 

a  constant-valued  random  vector.  Also,  let  h  map  R  into  R  and  denote 
by  D(h)  the  set  of  discontinuity  points  of  h. 
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(2.1.8)  LEMMA.  If  X ->X  and  Y  ->c,  than 

u  Q 

(X  .Y  )  ->  (X,c)  . 
n  n 

(2.1.9)  LEMMA.  If  X  —>X  and  either  h  la  continuous  or  P{XeD(h)}"0,  then 

n 

h(X  )  -o.  h(X)  . 

Q 

Lemma  (2.1.8)  [respectively  Lemma  (2.1.9)]  Is  a  special  case  of  Theorem  4.4 
[respectively  Theorem  S.l]  of  Billingsley  (1968). 

For  the  proof  of  Equation  (2.1.7),  we  first  note  that  the  standard 
c.l.t.  for  1.1. d. ,  mean  0,  finite  variance  random  variables  Implies  that 

l  n 

“777  Z  zk(f>  _>  M*1*  *  (2.1.10) 

an  '  k-1 

This  can  be  rewritten  as 

n1/2{?  (f )-r (f ) }/ [ (a/E{a, }) (E{o. }/a  ) ]  -o  N(0,1)  .  (2.1.11) 

u  x  in 

The  strong  law  of  large  numbers  (or  even  the  weak  law  for  that  matter) 
guarantees  that 

(E{a. }/a  )  ->  1  .  (2.1.12) 

i  n 

Lemma  (2.1.8)  applies  to  this  situation,  and  hence  if  we  let  Xn  denote 
the  left-hand  side  of  Equation  (2.1.11), 

(X  ,  (E{o.  }/<*_))  -»  (N(0,1)  ,1)  . 
n  in 

Mow  apply  Lemma  (2.1.9)  using  the  continuous  mapping  h(x,y)«x*y  to  conclude 


that 
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X  •(E{ot1}/cT  )  -»  N(0,1)  *1  .  (2.1.13) 

n  x  q 

Since  N(0,1)*1  has  Che  sane  distribution  as  N(0,1),  Equation  (2.1.13)  is 
the  sane  as  Equation  (2.1.7). 

Equation  (2.1.7)  provides  a  confidence  interval  for  r(f)  but  in 
general,  of  course,  the  "standard  deviation  constant"  a/ E{a^}  is  not 
available  and  nust  be  estimated.  The  most  straightforward  estimate  for 
o/E{a^}  is  sa/an,  end  the  strong  law  of  large  numbers  ensures  that 

(o/s  )  “>  1  . 
n 

Then,  by  the  same  argument  which  leads  to  Equation  (2.1.7),  we  obtain  from 
Equation  (2.1.10)  the  c.l.t. 

n1/2{?n(f)-r(f)}/(sn/an)  N(0,1)  . 

It  follows  that  for  0<y<l/2  the  interval 

Jn<f)  *  en(f)-'l-YV(v1/2)'  *n<f)+sl-yV(v1/2)]  1 

where  z1_^,*$~^(1-y)  and  $(•)  is  the  distribution  function  of  the 
standardized  normal  random  variable,  provides  an  asymptotically  valid 
100(l-2y)%  confidence  interval  for  r(f).  This  means  that 

lim  P{r(f)el  (f)}  -  l-2y  , 
w*00  n 

and  thus  when  simulating,  for  large  n  the  interval  ( f )  (having  random 

end  points)  surrounds  the  unknown  constant  r(f)  approximately  100(1-2y)Z 
of  the  time.  Other  point  and  interval  estimates  which  reduce  the  bias  of 


i  (£)  are  available;  see  Iglehart  (1975).  Note  that  the  Interval  I  (f)  - 
n  n 

Is  symmetric  about  the  point  estimate  £Q(f) •  and  that  the  half  length  of 

-1/2 

the  interval  is  n  times  a  multiple  (z^_^)  Che  estimate  of  the 
constant  a/E{ct^}.  Thus  as  n  increases,  the  length  of  the  Interval  converges 
to  0  and  the  midpoint  converges  to  the  true  value. 

The  procedure  just  described  for  obtaining  confidence  Intervals  when 
simulating  a  regenerative  process  is  fist  a  fixed  number  of  cycles.  For  a 
simulation  of  fixed  run  length  t,  the  procedure  (Crane  and  Iglehart 
(1975a) ,  p.  39)  is  the  same  except  that  statistics  are  computed  only  for 
the  (random  number  of)  cycles  completed  by  time  t.  Confidence  Intervals 
are  based  on  the  c.l.t.,  analogous  to  Eq.  (2.1.7) ,  that  as 

tl/2{en(t)(f)“r(f)}/to/(E{al})1/2]  M(0»15  » 

where  n(t)  is  the  number  of  cycles  completed  in  (0,t], 

Crane  and  Iglehart  (1975a)  have  also  shown  that  for  regenerative 
processes  possessing  more  than  one  sequence  of  regeneration  points,  with 
high  probability  the  resulting  confidence  Intervals  are  of  the  same  length, 
provided  that  the  length  of  the  simulation  run  is  large.  More  precisely, 
for  fixed  run  length  t,  if  I(t)  and  I*(t)  are  the  lengths  of  confidence 
Intervals  obtained  from  two  sequences  of  regeneration  points,  then 

lim  I(t)/I*(t)  -  1 
t-*» 


with  probability  one. 


Finally  ,  suppose  that  we  ara  lntarastad  in  estimating  r(f)  with  a 
100(l-2y)%  confidanca  interval  whose  half  length  is  1006%  of  r(f).  Than 
the  number  of  cycles  required  is  (approximately) 

n  -  (s1_Y/6)2ta/(r(f)E(a1})l2  . 

2 

The  first  factor  (z^_^/6)  is  independent  of  the  system  being  simulated. 

2 

From  the  second  factor  [a/(r(f)E{a^}> ]  it  is  apparent  that  some  systems 
are  inherently  more  dif fictile  to  simulate  than  others.  This  quantity 
provides  a  good  measure  of  the  length  of  simulation  required  for  a  fixed 
level  of  precision.  An  estimate  for  this  quantity  obtained  from  a  pilot 
run  can  be  used  to  determine  the  length  of  the  final  simulation  run. 

2.2.  Regenerative  Processes  in  Networks  of  Queues 

To  illustrate  the  ideas  of  Section  2.1,  consider  the  simple  closed 
network  of  queues  of  Figure  2.1.  A  fixed  number  of  jobs,  N,  circulate  in 
the  network  from  time  t*0.  Upon  completion  of  a  service  at  center  1,  a 
job  joins  the  tail  of  the  queue  in  center  2  for  6  service.  After  completion 
of  service  at  center  2,  the  job  joins  the  tall  of  the  queue  in  center  1. 
Neither  center  1  nor  center  2  service  is  subject  to  Interruption,  and  we 
assume  that  both  queues  are  served  according  to  a  first-come,  first-served 
(FCFS)  discipline.  We  complete  the  specification  of  this  network  of  queues 
by  making  the  following  probabilistic  assumptions: 

(i)  successive  a  service  times  form  a  sequence  of  1.1. d.  random 
variables  exponentially  distributed  (as  S^)  and  having  rate 
parameter  X^,  l.e.,  for  t20, 
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(il)  successive  8  service  times  form  a  sequence  of  l.l.d.  random 
variables  exponentially  distributed  (as  S^)  and  having  rate 
parameter 

(lii)  the  sequences  in  (i)  and  (li)  are  mutually  Independent. 

To  study  quantities  associated  with  the  lengths  of  the  queues  in 
centers  1  and  2,  for  tkO  let  X(t)  be  the  number  of  jobs  waiting  or  in 
service  at  center  1  at  time  t.  The  stochastic  process  X»{X(t):tkO}  has 
finite  state  space  E"{0,1, . . . ,N),  and  under  assumptions  (i),  (ii) ,  and 
(ill)  is  in  fact  a  continuous  time  Markov  chain.  It  is  easy  to  check  that 
the  Markov  chain  X  is  Irreducible  and  positive  recurrent.  Thus,  this 
"state  vector  process"  X  for  the  network  of  queues  is  a  regenerative 
process;  successive  returns  to  any  fixed  state  leE  are  regeneration  points 
for  the  process.  The  regenerative  process  £  has  a  "steady  state"  X,  and 
we  can  apply  the  regenerative  method  to  obtain  from  a  single  replication 
point  estimates  and  confidence  Intervals  for  quantities  of  the  form 
r(f)»E{f(X)}  for  some  function  f  having  domain  E.  Suppose  for  example, 
that  f  is  the  indicator  function,  1^  2  of  the  set  {l,2,...,N}. 

(The  indicator  function  1^  2  N} (x)  equals  1  if  xe{l,2, . . . ,N>  and 
equals  0  if  xl{l,2,. . . ,N},)  Then  r(f)"E{f(X)}  is  the  steady  state 
utilization  of  service  center  1,  i.e.,  the  (limiting)  probability  that 
center  1  is  busy.  If  f(x)«x,  then  r(f)  Is  the  steady  state  mean  number 
of  jobs  waiting  or  in  service  at  center  1. 

Mote  that  if  we  assume  that  one  of  the  service  time  random  variables 
has  an  exponential  distribution,  but  that  the  other  (say  S2)  has  a  finite 


■war*?  M  "  >  EJ  *"!Ji  .'TCg 
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mean  but  otherwise  arbitrary  distribution,  the  stochastic  process  X  is  no 
longer  a  continuous  tine  Markov  chain.  The  process,  however,  is  still 
regenerative;  the  successive  entrances  of  the  process  X  to  a  fixed  state 
lc{l,2, . . . ,N}  from  state  1-1  are  regeneration  points.  These  time  points 
correspond  to  completions  of  service  at  center  2  after  which  there  are  1 
jobs  at  center  1. 

Now  consider  the  network  of  queues  In  Figure  2.2,  formulated  (Lewis 
and  Shedler  (1971))  as  a  model  of  system  overhead  in  multlprogrammed 
computer  systems  operating  tinder  demand  paging.  The  Interpretation  of 
this  figure  differs  from  that  of  a  conventional  diagram  for  a  network  of 
queues  in  that  services  are  distinguished  from  the  servers  which  perform 
them.  In  Figure  2.2,  the  circles  represent  services  rather  than  servers. 
The  model  consists  of  two  sequential  stages,  the  a  stage  and  the  8  stage, 
in  a  loop.  Two  servers.  Interpreted  as  a  processor  and  a  data  transfer 
unit  (10  unit),  provide  service  to  a  constant  number,  N,  of  jobs 
(programs) ;  each  of  the  jobs  goes  through  both  stages  in  sequence  and  then 
returns  to  the  first  stage,  this  process  being  repeated  continuously. 

Within  the  a  stage,  a  job  receives  each  of  three  services  a^,  ,  and  <Xj, 

in  that  order  and  similarly,  within  the  8  stage,  a  job  receives  each  of 
three  services  8^,  Bj,  and  83,  1“  that  order.  A  service  can  be  provided 
only  by  the  10  unit,  and  each  of  the  other  services  can  be  provided  only 
by  the  processor.  We  assume  that  the  two  servers  can  provide  service 
concurrently,  subject  to  the  restriction  that  the  processor  cannot  provide 
a  8^  or  a  8^  service  while  the  10  unit  is  providing  a  82  service.  In 
addition,  we  assume  that  after  having  received  an  service,  a  job  moves 


waiting  station 


a-naga  0-stage 

(i)  Processor  provida*  ,  «j,  <*3, 0,  and  03  sarvicaa 

(ii)  I/O  unit  provida*  02  sarviea 

(iii>  No  j3, ,  03  sarviea  by  processor  during  0j  sarviea  by  I/O  unit 
<iv)  a, .  <*3. 0, .  0j.  03  sarvicas  not  intamiptabia 

(v>  atj  sarviea  hat  pra-amptiva  rasuma  typa  interruption  at  completion  of  02  service 

»  ' 

Figure  2.2.  System  overhead  model 
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Instantaneously  from  tha  a  stags  to  the  tall  of  the  queue  in  the  6  stage, 
and  after  having  received  a  service,  moves  Instantaneously  from  the  6 
stage  to  the  tall  of  the  queue  In  the  ct  stage. 

The  single  processor  provides  dp  a2,  Bp  and  6^  service  to  the 
several  jobs  in  the  netwrk.  Having  begun  an  oip  a^,  Bp  or  6^  service, 
the  processor  completes  the  service  without  interruption.  In  the  case  of 
the  service,  however,  Interruption  is  possible.  Interruption  of  an 
service  occurs  at  the  completion  of  a  concurrent  62  service.  After  some 
time,  the  a2  service  continues  from  the  point  In  the  service  at  which 
Interruption  occurred.  Thus,  the  "{^-complete"  interruption  of  an  a2 
service  is  of  the  preemptive-resume  type. 

At  the  completion  of  an  cip  a2,  a^.  Bp  or  &3  service  and  at  an 
Interruption  of  an  a2  service  (l.e.,  at  completion  of  $2  service),  the 
processor  chooses  the  next  service  to  be  provided  according  to  a  rule  of 
priority  as  follows  (see  Figure  2.3): 

(1)  if  there  is  a  job  waiting  for  B^  service,  begin  this  service; 

(11)  otherwise,  If  there  is  a  job  waiting  for  B^  service,  begin  this 
service  provided  that  B2  service  Is  not  in  progress; 

(111)  If  the  last  a  service  provided  was  a  completed  a2  service,  begin 
an  service; 

(iv)  If  the  last  a  service  provided  was  an  interrupted  a2  service, 
resume  the  a2  service; 

(v)  If  the  last  a  service  provided  was  an  service,  begin  an  a2 


service; 


28 


(vl)  if  the  last  a  service  provided  wee  an  service  and  if  che  queue 

in  the  a  stage  is  not  empty,  begin  an  service. 

If  no  claim  is  made  on  the  processor  according  to  the  rule  of  priority, 
it  remains  idle  until  the  completion  of  the  next  $2  service,  at  which  time 
che  rule  of  priority  is  Invoked  again.  We  assume  that  the  queue  in 
the  a  stage  and  the  queue  in  the  8  stage  are  both  served  according  to  a 
FCFS  discipline. 

To  complete  the  specification  of  the  model*  we  make  the  following 
probabilistic  assumptions : 

(1)  the  successive  ct^  [respectively  a^]  service  times  form  a  sequence 
of  1.1. d.  random  variables;  the  random  variable  [respectively 
a^]  has  a  finite  mean,  but  otherwise  arbitrary  distribution; 

(11)  the  successive  83  [respectively  82»  63]  service  times  form  a 
sequence  of  l.i.d.  random  variables;  the  random  variable  8^ 
[respectively  62*  63]  has  a  finite  mean,  but  otherwise  arbitrary 
distribution; 

(ill)  the  successive  a2  service  times  form  a  sequence  of  l.i.d. 
exponentially  distributed  random  variables; 

(iv)  the  sequences  in  (1),  (li)  and  (ill)  are  mutually  independent. 

Quantities  of  Interest  in  this  model  include  the  limiting  utilization 
of  the  10  unit  (i.e.,  the  long-run  expected  fraction  of  time  that  it  is 
busy  as  opposed  to  being  idle) ,  and  the  long-run  expected  fractions  of 
time  that  the  processor  provides  each  of  the  services  ot^,  c^,  ou ,  8^,  and 


S^i  Formally,  denote  the  total  amount  of  [reapaetively  a^,  a-y  g, , 

82,  831  sarvlca  that  haa  occurrad  In  tha  ayatam  during  tha  tlma  intarval 
(0,t]  by  A^t)  [raapectivaly  A2(t),  AjCt),  Bj^t),  B2(t),  B3(t)].  Than, 

In  tarma  of  thaaa  procaaaaa  of  cumulatad  aarvlca  timaa,  for  i»l,  2,  and 
3,  wo  wlah  to  aatlmata 

E{A. (t) } 

11m - t - 

and 

E{B . (t) } 

lim - * -  .  (2.2.1) 

provldad  that  tha  llmlta  azlat . 

Tha  daflnltlon  of  an  approprlata  ayatam  atata  for  a  simulation  of  thla 
modal  la  parhapa  aomawhat  laaa  apparont  than  In  tha  previous  example.  It 
la  dear,  however,  that  It  la  necasaary  to  take  into  account  more  than 
Juat  tha  number  of  jobs  waiting  or  In  aarvlca  at  one  of  the  stage8.  Tha 
additional  Information  that  wa  uaa  la  tha  kind  of  service  which  is  being 
provided  by  tha  processor  or  (according  to  tha  rule  of  priority)  Is  to  be 
provldad  next.  If  at  time  t  tha  processor  is  providing  service  (l£l£3) 
and  there  are  n  jobs  (Oin<fl)  waiting  or  In  aarvlca  In  the  8  stage,  we 
define  tha  system  state  X(t)  to  be  (n,i).  Otherwise,  at  time  t  the 
processor  is  Idle  or  la  providing  8^  or  83  service;  In  this  case  wa  define 
X(t)  to  be 

(1)  (n,0)  if  there  are  n  jobs  (OSnSN)  waiting  or  in  service  in  tha 

6  stage  and  tha  next  a  service  to  be  provldad  is  the  resumption 


of  an  a2  service; 


(11)  (n,l)  If  the  next  a  service  Is  an  service;  and 

(111)  (n,2)  If  the  next  a  service  Is  the  start  of  an  oig  service. 

The  force  of  this  state  definition  Is  that  the  stochastic  process 
g-{X(t) :t20)  Is  a  regenerative  process  In  continuous  time.  To  see  this, 
consider  the  increasing  sequence  of  time  points  { :k£0)  at  which  either: 

(1)  a  service  has  just  been  completed  and  the  served  job  has  moved 
to  the  a  stage  queue,  or 

(11)  after  a  time  point  described.  In  (1)  at  which  the  $  stage  queue  Is 
empty,  the  next  job  appears  at  the  8  stage  for  service. 

Mow  consider  the  subeequence  { t^}  of  the  {t^}  at  which  the  system  enters 
state  (1,1).  At  such  a  time  point,  there  is  one  job  In  the  B  stage  and 

the  next  a  service  to  be  provided  Is  an  service.  These  t£  are 

regeneration  points  for  the  process  g.  Mote  that  the  process  (X(tfc)  :lciO} , 
l.e.,  the  process  g  observed  at  the  epochs  {t^} ,  Is  an  Irreducible, 
aperiodic  finite  state  discrete  time  Markov  chain.  In  addition,  It  Is 
easy  to  check  from  the  basic  assumptions  of  the  model  that,  given  the 
state  of  the  system  at  time  points  t^  and  t^g,  the  time  Interval 
Is  a  random  variable  whose  distribution  does  not  depend  on  the  times 
between  previous  time  points  tg  (isk)  or  the  state  of  the  system  at  the 
time  points  tg  (i<k) ) . 

The  regenerative  structure  of  the  process  X  provides  a  basis  for 
estimation  of  the  quantities  defined  in  Equation  (2.2.1)  even  though  they 
cannot  be  expressed  In  the  form  E{f(X)}.  The  first  observation  Is  that 

the  sequences 
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{ ^ A± < c lc-4-1  > —  A±  ( C lc)  ’  tk+l”tk) :k21* 

and 

tkH_tk):k*1*  (2.2.2) 

consist  of  1.1. d.  random  vectors.  In  addition,  the  ratio  formulas 


E{A. (t) } 

lim - 1 - 

t-*» 


E{  Ai  ( t^i)  “Ai  ( } 


and 


E{B.(t)> 

11m - i - 

fK»  C 


E{Bl(tk+i>-Bl(ti)} 

Eltk+l-tk> 


(2.2.3) 


hold,  provided  that  all  relevant  expectations  are  finite.  The  ratio 
formulas  can  be  established,  for  example,  by  direct  application  of  a  basic 
limit  theorem  for  cumulative  processes  (Smith  (1958),  p.  263);  the  processes 
{A.^(t):tkO}  and  {B^(t):t*0}  are  cumulative  processes  with  respect  to  the 
regenerative  process  JC.  Given  Equations  (2.2.2)  and  (2.2.3),  the  arguments 
for  the  standard  regenerative  method  apply  and  yield  point  estimates  and 
confidence  Intervals  for  the  quantities  of  interest. 


This  somewhat  complex  model  Illustrates  the  following  remark  about 
network  of  queues  used  as  models  for  computer  and  communication  systems. 
For  quantities  associated  with  queue  lengths,  it  is  often  possible  to 
define  an  appropriate  state  vector  and  in  a  fairly  straightforward  manner, 
to  establish  the  existence  of  regeneration  points  and  the  applicability 
of  the  regenerative  method.  Frequently  this  is  possible  under  fairly 
general  distributional  assumptions.  The  key  to  showing  the  applicability 


of  the  regenerative  method,  as  in  the  example  above,  is  often  the  discovery 
of  a  familiar  stochastic  process  (e.g.,  discrete  time  Markov  chain) 
embedded  in  the  state  vector  process. 

There  are,  however,  examples  of  networks  of  queues  for  which,  even 
under  the  convenient  assumptions  of  distributions  with  a  Cox-phase 
(exponential  stage)  representation  for  service  and  interarrival  times,  it 
is  technically  quite  difficult  to  establish  the  applicability  of  the 
regenerative  method;  this  is  so  even  though  we  would  expect  the  required 
conditions  to  be  satisfied.  The  model  of  an  automated  tape  library 
proposed  by  Lavenberg  and  Slutz  (1975)  provides  one  such  example. 
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3.0.  CLOSED  NETWORKS  OF  QUEUES 

We  deal  with  closed  networks  of  queues  having  a  finite  number  of  jobs 
(customers),  N,  a  finite  number  of  service  centers,  s,  and  a  finite  number 
of  (mutually  exclusive)  job  classes,  c.  At  every  epoch  of  (continuous) 
time  each  job  is  in  exactly  one  job  class,  but  jobs  may  change  class  as 
they  traverse  the  network.  Upon  completion  of  service  at  center  1  a  job 
of  class  j  goes  to  center  k  and  changes  to  class  l  with  probability  p^ 
where  P^{p . .  ,  . :l£i,k£s,  l£j  ,JUc}  is  a  given  irreducible  Markov  matrix. 

Ij  *  JCA 

At  each  service  center  jobs  queue  and  receive  service  according  to  a  fixed 
priority  scheme  among  classes;  the  priority  scheme  may  differ  from  center 
to  center.  Within  a  class  at  a  center,  jobs  receive  service  according  to 
a  fixed  queue  service jdiscip line,  e.g.,  first-come ,  first-served  (FCFS). 
Note  that  in  accordance  with  the  matrix  P,  some  centers  may  never  see  jobs 
of  certain  classes.  Only  one  job  can  receive  service  at  a  center  at  a 
time;  i.e.,  the  centers  are  single  servers.  According  to  a  fixed  procedure 
for  each  center,  a  job  in  service  may  or  may  not  be  preempted  if  another 
job  of  higher  priority  joins  the  queue  at  the  center. 

3.1.  Probabilistic  Assumptions 

The  marked  job  method  discussed  in  Section  4  for  passage  time 
simulation  applies  to  networks  of  this  kind  in  which  all  service  times  in 
the  network  are  mutually  independent,  and  at  a  center  have  a  distribution 
with  a  Cox-phase  representation  (Cox  (1955)),  i.e.,  consisting  of  a 
sequence  of  exponential  stages;  see  Figure  3.1.  We  permit  parameters  of 
the  service  time  distribution  to  depend  on  the  service  center,  the  class 
of  job  being  served,  and  the  "state"  of  the  entire  network  as  defined  below. 
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Each  service  time  distribution  has  its  own  finite  nuinber  of  stages, 
say  n.  Realization  of  a  service  time  is  as  a  sum  of  a  random  number  (21) 
of  exponentially  distributed  times.  Within  the  jth  stage  (l£j£n-l),  an 
amount  of  service  exponentially  distributed  with  rate  parameter  accrues 
with  probability  1-b^  service  is  complete,  and  with  probability 
additional  service  exponentially  distributed  with  rate  parameter 
accrues.  The  density  function  of  the  resulting  service  time  has  rational 
Laplace  transform 


where  a^-1  and  for  l<j£n,  a^-b^. . .b^_^.  The  class  of  density  functions 
having  Laplace  transform  of  this  form  includes  hyperexponential  and 
mixtures  of  Erlang  densities;  see  the  Appendix  of  Gelenbe  and  Muntz  (1976) 
for  a  discussion.  Note  chat  we  exclude  the  case  of  zero  service  times 
occurring  with  positive  probability. 

In  connection  with  this  class  of  service  time  distributions,  it  is 

clear  that  if  b  ■b-"..."b  ."l  and  if  X,-X„-...«X  ,  then  the  resulting 

x  £  n-i  l  £  n 

service  time  has  an  Erlang-n  distribution,  l.e.,  a  gamma  distribution  with 
integral  shape  parameter.  It  is  less  obvious  but  can  be  shown  (by 
considering  Laplace  transforms  of  the  density  functions)  that  if 
and 


bl  -  1  -  pl  -  P2^2/Xl  * 

then  the  Cox-phase  form  is  a  representation  for  the  hyperexponential 
density 
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-A.t  -A,t 

p^A^*  +  ?2^2*  * 

concentrated  on  t*0,  where  end  p2  ere  nonnegetlve  end  Pj+p2*l.  The 
corresponding  result  for  three  stages  is  that  if  *i>^2>X3* 

bx  -  1  -  Px  -  (PjAj/A^  -  (p3A3/Ax)  , 

and 

b2  -  1  -  IP2+{p3A3(A1-A3)}/{A2(A1-A2}l/[p2+{p3(A1-A3)/(A1-A2)}]  , 

then  the  Cox-phase  form  is  a  representation  for  the  hyperexponential 
density 


-At  -A2t  -At 

plXle  +  P2X2®  +  P3X3*  * 

where  p^,  p2,  and  p3  are  nonnegative  and  p^+p2+p3«l. 

3.2.  State  Vector  Definition 

In  order  to  characterize  the  state  of  the  network  at  time  t,  we  let 

S^(t)  denote  the  class  of  the  job  receiving  service  at  center  i  at  time 

t,  where  1*1,2 . .  by  convention  S^tJ^O  if  at  time  t  there  are  no  jobs 

at  center  i.  The  classes  of  jobs  serviced  at  center  i  ordered  by 

decreasing  priority  are  j2(i)  ,j2(i) , . . .  .j^^  (i)  ,  elements  of  the  set 

{l,2,...,c}.  Let  cf*\t) . (t)  denote  the  number  of  jobs  in  queue 

J1  Jk(i) 

at  time  t  of  the  various  classes  of  jobs  serviced  at  center  i,  1*1,2 . s. 

For  queue  lengths  of  jobs  of  various  classes  at  the  several  centers,  these 
state  variables  (together  with  the  stages-of-service)  would  suffice.  To 
deal  with  general  characteristics  of  passage  times,  however,  these  state 
variables  are  Inadequate.  An  apperently  minimal  state  vector  augmentation 
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Is  based  on  the  concept  of  a  marked  job.  The  idea  is  to  keep  track  of 
the  position  of  an  arbitrarily  chosen  marked  job  in  the  network  and  to 
measure  its  passage  times  during  the  simulation.  It  is  convenient  to 
think  of  the  N  jobs  as  being  completely  ordered  in  a  linear  stack  according 
to  the  following  scheme.  For  t*0,  we  define  the  vector  Z(t)  by 

Z(t)  -  (C.(1)  (t),C.(1)  (t) . cJ1)(t),S.(t);...i 

Jk(l)  Jk(l)-1  31 

C.(s)  (t),C.(s)  (t) . C.(s)(t),S  ft»  .  (3.2.1) 

Jk(s)  Jk(s)-1  J1 

The  linear  job  stack  then  corresponds  to  the  order  of  components  in  the 
vector  Z(t)  after  ignoring  any  zero  components.  Within  a  class  at  a 
particular  service  center,  jobs  waiting  appear  in  the  job  stack  in  FCFS 
order,  i.e.,  in  order  of  their  arrival  at  the  center,  the  latest  to  arrive 
being  closest  to  the  top  of  the  stack.  We  denote  by  N(t)  the  position 
(from  the  top)  of  the  marked  job  in  this  job  stack  at  time  t. 

We  associate  a  stage  of  service  with  each  job  in  the  network  as 
follows.  A  job  in  service  is  in  a  particular  stage  of  its  service  time 
distribution  at  that  center;  for  such  a  job,  this  is  the  associated  stage. 
For  a  job  in  queue,  the  associated  stage  is  the  stage  of  its  service  time 
distribution  that  is  to  be  provided  when  the  job  next  receives  service; 
typically  this  is  the  first  stage  of  service,  but  may  be  a  subsequent 
stage  if  the  job  has  been  preempted.  For  t£0,  we  define  the  vector  U(t) 
by 

U(t)  -  (U1(t),...,UN(t))  ,  (3.2.2) 
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where  (t)  is  the  stage  of  service  eesocieted  with  the  jth  job  in  the 
linear  stack  of  jobs,  and  for  tiO  take  as  the  state  vector  of  the  network 
of  queues  the  vector 


X(t)  -  (Z(t),N(t),U(t))  .  (3.2.3) 

For  any  service  center  1  that  sees  only  one  class  of  job,  l.e.,  such  that 

k(i)-l,  It  is  possible  to  simplify  the  state  vector  by  replacing 

cf1)  (t),S.(t)  by  Q. (t),  the  total  number  of  jobs  at  center  1.  Mote  that 
Jk(i)  1  1 

the  state  vector  definition  does  not  explicitly  take  into  account  that 
the  total  number  of  jobs  in  the  network  is  fixed.  In  the  case  of  complex 
networks,  the  use  of  this  resulting  somewhat  larger  state  space  facilitates 
generation  of  the  state  vector  process;  for  relatively  simple  networks, 
it  may  be  desirable  to  remove  the  redundancy. 

3.3.  Definition  of  Passage  Times 

Given  a  particular  closed  network  of  queues,  we  must  specify  the 
passage  time  of  interest.  This  can  be  done  in  terms  of  the  arbitrarily 
chosen  marked  job  of  Section  3.2,  by  means  of  four  subsets  (A^  A^,  B^,  Bg) 
of  the  state  space,  E,  of  the  stochastic  process  £-{X(t) :t*0}.  The  sets 
A^,  A2  [respectively  B^,  B2]  jointly  define  the  random  times  at  which 
passage  times  for  the  marked  job  start  [respectively  terminate] .  The  sets 
A^,  A2>  Bx  and  B2  in  effect  determine  when  to  start  and  stop  the  clock 
measuring  a  particular  passage  time  of  the  marked  job. 

It  is  convenient  to  introduce  the  jump  times  {xQ:n20}  of  the  process  £* 
We  set  Tq-0  and  have  tq<t1<. . .  with  probability  one.  Since  the  state 
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space  E  of  ^  la  finite,  chare  can  be  no  finite  accumulation  points  for 
the  Tq  and  t  -**.  For  k,nkl,  we  require  that  the  sets  A^,  A2,  and  B2 
satisfy : 

if  X^^sA^  X(fn)€A2,  X(T  n-l+k^  cAl  and  X^Tn+k^£A2* 
then  “d  x^im)£B2  for  some  0<m*k; 

and 

if  X(tn_1)«B1,  X(Tn)«B2,  X(xq_ 

then  X(T  _1+m)eAi  and  X(T^)m)eA2  for  some  0£m<k  .  (3.3.1) 

These  conditions  ensure  that  the  start  and  termination  times  for  the 
specified  passage  time  strictly  alternate. 

In  terms  of  the  secs  A^  Aj,  B^  and  B2,  we  define  two  sequences  of 
random  times,  {S^ :jaO}  and  {T^ :jil},  where  Sj_x  is  the  start  time  of  the 
jth  passage  time  for  the  marked  job  and  is  the  termination  time  of  this 
jth  passage  time.  Assuming  that  the  initial  state  of  the  process  X  is  such 
that  a  passage  time  for  the  marked  job  atazts  at  t“Q ,  let 

S0"° 

Sj  -  inf  :X(x^)eA2,  X(xk)eAx  for  some 

and  k<n} ,  jil  , 
and 

Tj  -  inf{T0>Sj_1:X(Tn)eB2,  X(xk)eB1  for  some 
x^iS^  and  k<n},  jil  . 


(3.3.2) 
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Then  Che  jth  passage  tins  for  the  marked  job  is  simply  J2l. 

Quite  often  (and  in  particular  in  the  two  examples  discussed  in  Section  5) , 
the  value  of  k  in  the  definitions  of  Sj  and  can  be  set  equal  to  n-1, 
a  considerable  simplification.  In  ganaral,  however,  we  allow  k<n.  For 
response  times  (complete  circuits),  A^-B^  and  consequently 

for  all  jkl. 

The  following  example  illustrates  specification  of  passage  times  in 
terms  of  the  sets  A^,  ,  B^,  and  B^. 

(3.3.3)  EXAMPLE.  Consider  the  closed  network  of  Figure  3.2  having 
three  service  centers.  Upon  completion  of  a  [respectively  B]  service  at 
center  1  [respectively  center  2],  each  of  the  N  jobs  joins  the  tail  of  the 
queue  in  center  3.  In  accordance  with  a  binary-valued  variable  ty,  after 
completion  of  y  service  at  center  3,  a  job  joins  the  tall  of  the  queue  in 
center  1  (when  iji-l)  or  joins  the  tail  of  the  queue  in  center  2  (when  i|>"0). 
Assume  that  each  of  the  queues  is  served  according  to  a  FCFS  discipline. 

Also  assume  that  routing  variable  is  a  random  variable  independent  of 
service  times  in  the  network,  and  that  at  each  of  the  centers,  service 
times  are  l.l.d.  exponential  random  variables. 

Suppose  that  the  passage  time  of  Interest  starts  when  a  job  joins  the 
center  1  queue  (upon  completion  of  service  at  center  3)  and  terminates  when 
the  Job  joins  the  queue  in  center  3.  Here  there  is  effectively  a  single 
job  class.  Denoting  by  p  the  (independent)  probability  that  the  binary 
routing  variable  <|>  takes  the  value  1,  in  the  routing  matrix  P_,  p. .  .«p9 

"  IX I XX) JX 
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p31,ll“p*  P31,2l“1~P’  *nd  #11  oth*r  pij ,klm°‘  For  ti0  and  i_1>2»  ^  3» 
we  define  Q^(t)  to  be  the  number  of  job*  waiting  or  in  service  at  center  i. 


X(t)  -  (Z(t) ,N(t))  , 


where  Z(t)*(Q^(t) ,Q2(t) ,Q3(t))  and  N(t)  is  the  position  of  the  marked  job 
in  the  associated  job  stack.  The  state  space  E  of  the  process  £  is 


E  -  {(i,j ,k,n) :0Si,j,ksN;  i+j+k-N;  lSoSN}  . 


For  this  passage  time 


A1  -  {(i,j ,k,n)€E;k£l;  n-N}  , 
A2  *  {(i,j ,k,n)eE:i2l;  n-l}  , 
-  {(i,j  ,k,n)€E:iil;  n«i}  , 


B2  -  {(i,j,k,n)«E:k2l;  n-i+j+1} 


How  assume  that  the  queues  in  center  1  and  center  2  are  served 
according  to  a  FCFS  discipline,  but  that  at  center  3,  Jobs  that  join  the 
queue  after  completion  of  service  at  center  1  have  priority  over  those  that 
join  the  queue  after  completion  of  service  at  center  2.  Suppose  that  the 
passage  time  (indicated  by  F  in  Figure  3.2)  of  interest  starts  when  a  job 
joins  the  queue  itr  center  1  and  terminates  at  the  completion  of  service 
to  the  job  at  center  3. 


To  deal  with  this  passage  time,  we  use  two  job  classes:  center  1  sees 
only  jobs  of  class  1,  center  2  sees  only  Jobs  of  class  2,  and  center  3  sees 


43 


jobs  of  class  1  and  class  2.  In  Che  routing  matrix  P,  p^  21^22  32*^' 
p31, 11*^32, 11*^'  p 31, 22*^32, 22*1'"'  “d  *U  oth*r  flj.kl*0'  Fot  «“>• 

Z(c>  -  (Q1(t),Q2(t),c”)(t),cP>(t).S3(t))  , 

(3)  (2) 

where  7 (t)  [respectively  (t)J  is  the  number  of  class  1  [respectively 
class  2]  jobs  in  queue  at  center  3  at  time  t,  and  S^(t)  is  the  class  of  job 
in  service  at  center  3  at  time  t.  (S^(t)«0  if  at  time  t  there  are  no  jobs 
at  center  3).  The  state  space  E  of  the  process  X*{X(t);t^O}  is 

E  -  {(i,j,0,0,0,n):0si,jsN;  i+j-N;  lSnSN}  u 

{ (i,j  ,lc,A,m,n)  :0£i,j  ,k,l£N;  i+j+k+E-N-1;  lSm£2;  lswtfl} 

For  the  passage  time  P, 

^  -  {(i,j ,k,4,m,n)«E:m-l;  n-N}  , 

A2  -  {(i,j ,k,A,m,n)«E:ikl;  n-l}  , 

B1  *  A1  * 

and 

B2  -  A2  °  {(i,j ,k,A,m,n)eE;j*l;  n-i+l}  . 
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4.0.  THE  MASKED  JOB* METHOD 

We  base  the  estimation  of  general  characteristics  of  passage  times 
in  a  closed  network  of  queues  of  Section  3  on  the  measurement 
of  passage  times  for  a  typical  job,  the  marked  job  discussed  above.  It 
is  Intuitively  clear  and  is  shown  in  Appendix  1  that  the  sequence  of  passage 
times  for  any  other  job  (as  well  as  the  sequence  of  passage  times, 
irrespective  of  job  identity,  in  order  of  start  or  termination)  converges 
in  distribution  to  the  same  random  variable  as  the  sequence  of  passage 
times  for  the  marked  job.  It  follows  that  we  can  estimate  general 
characteristics  of  passage  times  in  the  network  by  simulation  of  the 
process  X*{X(t) :t20}  defined  by  Equations  (3.2.1),  (3.2.2)  and  (3.2.3). 

We  shall  see  that  it  is  possible  to  use  the  regenerative  method  (applied 
to  an  appropriate  discrete  time  regenerative  process)  to  obtain  strongly 
consistent  point  estimates  and  asymptotically  valid  confidence  Intervals 
for  passage  time  characteristics. 

We  denote  by  Xq,  n20,  the  state  of  the  system  when  the  (n+l)st  passage 
time  of  the  marked  job  starts.  For  jsl,  let  Pj  be  the  jth  passage  time 
for  the  marked  job  and  take  the  quantity  of  interest  in  the  simulation  to 
be 


r(f)  -  E{f (X,P) }  ,  (4.0.1) 

where  f  is  a  real-valued  (measurable)  function  and  (X  ,P  ,.)“>(X,P).  For 

n  nTX 

example,  to  estimate  E{P),  we  take  f(x,p)ap;  to  estimate  P{Pst},  we  take 


is  the  indicator  function  of  the  set  [0,t]. 
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4.1.  Simulation  for  Passage  Times 

Based  on  the  measurement  of  passage  times  for  the  marked  job,  we  obtain 
estimates  of  the  quantity  r(f)  as  follows. 


(4.1.1)  ALGORITHM.  (Marked  Job  Method) 

1.  To  serve  as  a  return  state,  select  a  star*  of  the  system,  1q,  at 
which  a  passage  time  of  the  marked  job  starts.  Begin  the 


2. 


3. 


4. 


simulation  with  X(0)-1q. 

Carry  out  the  simulation  of  %  for  a  fixed  number,  n,  of  cycles 
(having  random  length)  defined  by  the  successive  returns  to  the 
state  1q. 

In  each  cycle  measure  all  passage  times  of  the  marked  job.  Denote 
the  number  of  passage  times  of  the  marked  job  observed  in  the 
kth  cycle  by  a^(kkl)  and  compute  the  quantity 


V1 


Vf> 


n-g 


f(Xn’Vl>  * 


(4.1.2) 


k-1 


where  gg-0  and  g^-o^+. .  ,+ct^. 

Take  as  a  point  estimate  (based  on  n  cycles)  of  r(f)  defined  by 
Equation  (4.0.1)  the  quantity 


5. 


where 


and 


?  (f)  -  Y  (f) /a  , 
n  n  n 


Y  (f)  -  n"1  E  Y.  (f) 
*  k-1  k 


-1 


(4.1.3) 


aQ  -  ° "  £  \  • 
n  k-l  * 

Take  as  a  100(l-2y)Z  confidence  Interval  (based  on  n  cycles)  for 
r(f)  the  Interval 


ai 
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Vf> 


[Vt)-‘lVa/(vl/J)' 


(4.1.4) 


Here  z^  ^(l-y) ,  where  $(*)  Is  the  distribution  function  of 
e  standardized  (mean  zero,  variance  one)  normal  random  variable, 
and  s&  is  the  quantity 


where 


and 


*„  •  !*n-2Vf>*i2*(V(>)2»22}1/2 

s..  -  (n-l)-1  £)  (Y.  (f)-Y  (f))2  , 
k-1 

*22  ■  t  <«A)2 

*12  -  fa-u'1 1  • 


Justification  for  this  passage  time  simulation  method  appears  in  the  next 
section. 


4.2.  The  Underlying  Stochastic  Structure 

With  knowledge  of  the  matrix  P  and  the  parameters  of  the  exponential 
stage  service  times,  we  can  carry  out  the  simulation  (i.e.,  generation  of 
the  process  X»{X(t) :t£0)  as  defined  in  Section  3.2)  in  a  straightforward 
manner.  Note  that  when  all  service  times  are  exponentially  distributed, 
the  vector  U(t)  in  the  state  vector  given  by  Equation  (3.2.3)  can  be 
omitted,  resulting  in  a  smaller,  less  complex  state  space. 

The  force  of  the  assumptions  of  Cox-phase  service  times  and  Markovian 
routing  in  the  closed  networks  of  queues  of  Section  3.1  is  contained  in 
Theorem  (4.2.1). 
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(4.2.1)  THEOREM.  The  stochastic  process  X-{X(t)  :tkQ}  defined  by 
Equations  (3.2.1),  (3.2.2)  and  (3.2.3)  is  a  continuous  time  Markov  chain 
with  a  finite  state  space,  E. 


The  method  of  Section  4.1  for  estimation  of  general  characteristics 
of  passage  times  relies  on  the  measurement  of  passage  times  for  a  typical 
job,  the  marked  job  as  introduced  in  Section  3.2.  As  in  Section  3.3, 
specification  of  the  passage  time  of  interest  is  in  terms  of  the  marked 
job,  by  means  of  the  four  subsets  (A^,  A^,  B^,  B^)  of  E.  We  define  the 
two  sequences  of  random  times,  {S^:j20}  and  (T^  :  j2l) ,  where  is  the 

start  time  of  the  jth  passage  time  for  the  marked  job  and  is  the 
termination  time  of  this  jth  passage  time.  For  an  initial  state  of  the 
Markov  chain  X  such  that  ,a  passage  time  for  the  marked  job  begins  at  t**0. 


s0.° 

S.  -  inf{x  2T  :X(T  )cA«,  X(T.)eA,  for  some 
j  n  j  n  l  k  1 


and  k<n},  J21 


T.  »  inf{x  >S.  ,:X(x  )eB_,  X(x,  )eB.  for  some 
j  nj-i  n  i  kl 


x^aS^^  and  k<n},  jsl  , 


(4.2.2) 


where  {xQ:n20}  are  the  jump  times  of  the  Markov  chain  X.  Then  the  jth 
passage  time  for  the  marked  job  is  Pj*Tj-Sj_^,  j2l. 


Let  X  denote  the  state  of  the  continuous  time  Markov  chain 
n 

X-{X(t):t20}  when  the  (rrt-l)st  passage  time  of  the  marked  job  starts: 

X  -X(S  ),  nkO.  Since  X  is  a  Markov  chain  and  {S  :nkO}  are  stopping  times 
n  n  ~  n 
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for  the  chain,  {X  :niO}  is  a  discrete  time  Markov  chain  with  finite  state 

Q 

space  Aj.  For  a  discussion  of  stopping  times,  see  finlar  (1975a),  p.  239. 

Throughout,  we  shall  assume  that  the  process  {X  :niO}  is  irreducible  and 

n 

aperiodic.  Furthermore,  it  is  easy  to  check  that  the  process  { (Xq,Sq) :n^0) 
satisfies 

F*Xn+l"J'  WS«StlX0 . V  S0 . Sn^ 

•  PtVl-J  •  WS»StlXa> 

with  probability  one  for  all  nkO,  j  ,  and  tkO. 

(4.2.3)  THEOREM.  The  stochastic  process  {(X  ,S  ) :n£0}  is  a 

n  n 

Markov  renewal  process. 

This  follows  directly  from  the  definition  of  a  Markov  renewal  process; 
see  Cinlar  (1975b)  or  Cinlar  (1975a),  p.  313.  The  basic  data  for  this 
Markov  renewal  process  is  the  semi-Markov  kernel  over  A^ , 

K-{K(i,j;t) :i,JeA2,  tkO},  defined  by 

KCl.iJt)  -  Vrsn“lvl}  • 

While  the  kernel  K  is  normally  given  in  the  analysis  of  a  Markov  renewal 

process,  for  the  network  of  queues  passage  time  problem,  K  is  virtually 

impossible  to  calculate.  Thus  from  this  point  of  view,  the  only  hope  for 

studying  the  passage  time  in  question  is  to  generate  sample  paths  of 

{(X  , S  ) :nkO)  via  simulation  of  the  Markov  chain  X. 
n  n 

As  in  Section  4.0,  let  f  be  a  real-valued  (measurable)  function  with 
domain  A2*R+,  where  R^tO,80),  and  define  the  quantity  r(f)  according  to 
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r (f)  -  E{f(X,P)}  .  (4.2.4) 

Mow  we  select  a  fixed  state,  ig,  from  A begin  the  simulation  of  2  with 
X(0)-i0,  and  carry  out  the  simulation  of  X  in  cycles  defined  by  the 
successive  returns  to  state  ig.  Let  denote  the  length  (in  discrete 
time  units)  of  the  kth  cycle  of  {XQ:niO}  and  define  8g-0,  and  .  •+ot^» 

k2l. 


Theorem  (4.2.5)  follows  from  Theorem  (4.2.3)  and  the  definition  of  a 
regenerative  process.  It  comprises  the  key  observation  which  leads  to 
point  and  interval  estimates  for  r(f). 

(4.2.5)  THEOREM.  The  stochastic  process  {(X  ,P  ,, ) :n*0}  is  a  regenerative 

n  nrl 

process  in  discrete  time  with  regeneration  points  (B^rk^O). 

Note  that  the  regeneration  points  8^  are  not  the  times  of  return  to  a 

fixed  state  of  the  process  {(X  ,P  .,) :naO}. 

n  ttri 

For  k^l ,  we  now  define 

V1 

Yk(f)  -  E  f(Xn,Pn+1)  .  (4.2.6) 

n  6k-l 

Since  the  8^  are  regeneration  points  for  {(X^.P^^) :n£0),  we  have  (cf. 

Crane  and  Iglehart  (1975),  Proposition  2.1)  the  following  result. 

(4.2.7)  THEOREM.  The  pairs  of  random  variables  {(Yk(f) ,0^) :k2l}  are 
independent  and  identically  distributed. 

The  regenerative  property  guarantees  (Miller  (1972))  that  as  n-*00 

<VW  *  «•»>  > 
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where  ->  denotes  convergence  in  distribution,  i.e.,  there  ere  rendoa 
variables  X  and  P  such  that  as  n-*=°,  the  joint  distribution  of  (X  ,P^^) 
converges  to  the  distribution  of  (X,P).  For  the  function  f  (appearing  in 
Equation  (4.2.4))  in  the  definition  of  r(f),  let  Df  denote  the  set  of 
discontinuities  of  f.  Assuming  that  P{ (X,P)«Df}«0  and  using  Lemma  (2.1.9) 
it  follows  that 

f(Xn,Pttfi)  ->  f(X,P)  (4.2.8) 

as  rr *».  The  final  step  is  to  establish  a  ratio  formula  for  E{f(X,P)} 
which  makes  it  possible  to  apply  the  regenerative  method  to  estimation  of 
passage  times;  this  follows  from  the  general  result  for  regenerative 
processes  (cf.  Crane  and  Iglehart  (1975),  Proposition  A. 3).  A  direct 
proof  not  requiring  the  key  renewal  theorem  is  in  Appendix  2. 

(4.2.9)  THEOREM.  Assume  that  E{ | f (X,P) | }<®.  Then 

E{f(X,P)}-E{Y1(f)}/E{a1>  , 

where  Y^f)  is  given  by  Equation  (4.2.6). 

With  the  ratio  formula  of  Theorem  (4.2.9)  and  the  fact  that  the  pairs 

of  random  variables  { (Y^(f) .o^) :kil}  are  independent  and  identically 

distributed  (Theorem  (4.2.7)),  we  can  apply  the  regenerative  method  to 

{ (X  ,P  ):n20}  to  obtain  point  and  Interval  estimates  for  r(f). 
n  n+i 


5.0.  EXAMPLES  AND  SIMULATION  RESULTS 


To  illustrate  the  marked  job  method  of  the  previous  sections  for 
estimation  of  passage  times  in  closed  networks  of  queues,  we  consider  two 
examples.  Descriptions  of  the  networks  and  state  vector  definitions  appear 
in  this  section,  along  with  illustrative  numerical  results. 

5.1.  A  Closed  Network  of  Queues 

The  first  example  is  relatively  simple;  see  Figure  5.1.  Upon  completion 
of  a  service  at  center  1,  in  accordance  with  a  binary-valued  variable  ty, 
the  job  joins  the  tall  of  the  queue  in  center  1  (when  ij^l)  or  joins  the 
tail  of  the  queue  in  center  2  for  $  service  (when  After  completion 

of  6  service  at  center  2,  the  job  joins  the  tail  of  the  queue  in  center  1. 

Neither  center  1  nor  center  2  service  is  subject  to  interruption.  We 
assume  that  both  queues  are  served  according  to  an  FCFS  discipline.  The 
limiting  response  time  of  interest  (denoted  by  R)  for  a  job  is  the  time 
measured  from  entrance  into  the  center  1  queue  upon  completion  of  a  center 
2  service  until  the  job  next  joins  the  center  1  queue.  Also  of  possible 
interest  in  this  model  is  the  limiting  passage  time  (denoted  by  P) 
corresponding  to  the  time  measured  from  entrance  into  the  center  1  queue 
upon  completion  of  i  center  2  service  until  the  job  next  joins  the  center 
2  queue. 

We  consider  passage  time  simulation  of  this  network  of  queues  under 
the  following  probabilistic  assumptions  t 

(i)  successive  a  service  times  form  a  sequence  of  i.i.d.  random 

I 

f 

variables,  exponentially  distributed  with  rate  parameter  y^;  j 

(il)  successive  0  service  times  form  a  sequence  of  i.i.d.  random 
variables,  exponentially  distributed  with  rate  parameter  yx ; 

1 


(i)  a  and  0  services  are  not  interruptable 

(ii)  Routing  determined  by  binary  valued  random  variable  ^ 

Figure  5.1.  Closed  network  of  queues 
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(Hi)  ip  is  a  Bernoulli  random  varlabla  and  valuas  of  ij)  at  auceaaaiva 

a  sarvica  completions  form  a  sequence  of  i.i.d.  random  varlablea; 
(iv)  tha  sequences  In  (i) ,  (ii)  and  (lii)  are  mutually  independent. 

In  this  network,  there  are  two  classes  of  jobs:  class  1  jobs  at 
center  1  and  class  2  jobs  at  center  2.  Since  each  center  sees  only  one 
class  of  job,  the  process  {Z(t):t20}  can  be  defined  by 

Z(t)  -  (Q1(t),Q2(t>) 

where 

Q^Ct)  •  number  of  jobs  waiting  or  in  service  at  center  1  at  time  t, 
and 

Q2(t)  ■  number  of  jobs  waiting  or  in  service  at  center  2  at  time  t. 


Now  let  N(t)  denote  the  position  of  the  marked  job  in  the  linear  job 
stack  corresponding  to  the  order  of  the  nonzero  components  of  Z(t).  Taking 
into  account  the  fixed  number  of  jobs  in  the  network  and  that  the  service 
times  are  exponentially  distributed,  the  process  X*{X(t) : t£0) ,  where 

X(t)  -  (Q^t)  •■(«))  , 

is  the  underlying  continuous  time  Markov  chain. 

For  this  modal  the  etate  spaca  E  of  X  is 

E  -  {(i,j)  :0HSN;  lSjsN}  , 

where  N  is  the  number  of  jobs  in  the  network.  In  the  special  case  of 
N«2  jobs,  Figura  5.2  shows  tha  possible  state  transitions  among  the 
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| 

N* (N+l)-6  states.  State  (0,1)  corresponds  to  both  jobs  in  canter  2  with 
the  narked  job  being  the  one  in  queue;  state  (1,1)  corresponds  to  one  job 
in  service  at  center  1  and  one  in  service  at  center  2,  the  marked  job 
being  the  former,  etc. 

The  Markov  chain  governing  change  of  job  class  at  service  completions 
can  be  taken  to  have  state  space  {1,2}  and  transition  matrix 

P  1-P 

l  - 

1  0 

where  we  denote  by  p  the  probability  that  the  binary  branching  variable  ij; 
takes  the  value  1. 

The  sets  and  (of  Section  4.2)  defining  the  start  of  response 
times  for  the  marked  job  are 

^  -  {(i,j)«E:j-N,  i<N} 

and 

^  -  {(i.j)eEsj-l,  i>0)  . 

Since  R  is  a  response  time,  Tj*S^  for  j2l,  B^-A^,  and  B^"^’  For  N-2  jobs, 

Aj-{(0,2),  (1,2)}  and  A2«{(1,1),  (2,1)};  see  Figure  5.2.  A  sample  path  of 
the  Markov  renewal  process  {(XQ,Sn) :nk0}  is  in  Figure  5.3.  For  the  passage 
time  P,  A1-{(0,2),  (1,2)},  A2«{(1,1),  (2,1)},  Bj-Ul.l),  (2,2)}  and 
B2-{(0,1),  (1,2)}. 

5.2.  A  Computer  System  Model 

The  second  example  is  mors  complex,  and  is  essentially  the  network  of 
queues  defined  by  Lavenberg  and  Shadier  (1976)  as  a  model  of  resource 


contention  In  the  so-called  "DL/I  component"  of  an  IMS  (Information 
Management  System)  data  base  management  computer  system;  see  Figure  5.4. 

The  interpretation  of  this  figure  differs  from  that  of  a  conventional 
queueing  network  diagram  in  that  services  are  distinguished  from  the 
servers  which  perform  them.  Thus,  circles  in  Figure  5.4  represent  services 
rather  than  servers.  Five  types  of  service,  denoted  by  otg,  a,,  a,,  au 
and  3  are  represented  in  this  model.  The  a  services  (a^,  a^,  and  a^) 
are  performed  by  a  single  server,  interpreted  as  the  processor,  and  the 
8  service  is  performed  by  a  single  server.  Interpreted  as  an  input-output 
(I/O)  unit.  We  assume  in  this  model  that  no  two  a  services  can  be  in 
progress  concurrently,  but  that  any  a  service  can  be  in  progress 
concurrently  with  a  3  service.  We  also  assume  that  each  of  the  a^,  a1 , 
a and  3  services  is  noninterrup table.  The  service  is,  however, 
lnterruptable  at  the  completion  of  a  3  service,  this  interruption  being 
of  the  preemptive-resume  type. 

A  fixed  number  of  jobs  circulate  in  the  network  from  time  zero.  At 
any  subsequent  instant  of  time,  a  job  either  is  receiving  service  or 
waiting  for  service  in  one  of  the  queues.  We  assume  that  each  of  the 
queues  is  served  according  to  an  FCFS  discipline.  Note  that  there  are 
two  queues  (denoted  by  q^  ^  and  q^  2)  f or  service  and  two  queues  (q^  ^ 

and  q^  2)  f°r  <*2  sarvice*  The  arrows  in  Figure  5.4  indicate  the  flow  of 
jobs.  There  are  two  branches  leaving  the  service  and  two  branches 
leaving  the  a ^  service;  these  branches  are  labeled  by  binary-valued 
variables  and  Upon  completion  of  an  or  service,  the  job  just 


(i)  Processor  renders  aQ,  a, .  a2  and  a3  services 
(iil  I/O  unit  renders  *3  service 

(iii)  aQ,  a, ,  a2  and  (3  services  are  not  interruptible 

(iv)  a3  service  has  pre-emptive  resume  type  interruption  at  completion  of  J  service 
(vl  Processor  scheduled  according  to  priority  ordering  of  queues  q0,  q1t ,  q,  2,  q2  , 

q2Jandq3 

<vi)  Pouting  determined  by  binary  valued  random  variables  J/ ,  and  ^2 


Figure  5.4.  OL/I  component  model 
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served  follows  the  brench  with  e  lebel  having  value  1.  Thus,  e.g.,  If 
upon  completion  of  an  service  the  current  value  of  Is  zero,  the  job 
just  served  joins  queue 


An  epoch  of  completion  of  any  a  service,  or  an  epoch  of  completion  of 
a  6  service  at  which  either  no  a  service  Is  in  progress  or  an  service 
is  in  progress.  Is  called  a  scheduling  decision  epoch.  At  such  a  time 
point ,  the  next  processor  service  is  determined  by  a  processor  scheduling 
algorithm.  Upon  completion  of  a  service,  the  served  job  immediately  joins 
the  next  queue  on  its  route.  The  next  processor  service  is  the  service 
having  highest  priority,  this  priority  being  determined  by  a  total  ordering 
of  queues  qQ,  q^^  ^ ,  q^^  2>  92  1*  q2  2  ^  q3*  The  Processor  scheduling 
algorithm  employs  the  total  ordering  q^^  2,  q2  ^  q^,  q2,2*  q0’  q3 
(highest  Co  lowest  priority )« 


The  limiting  passage  time  of  interest  in  the  DL/I  component  model  Is 
indicated  by  R  in  Figure  5.4.  It  is  interpreted  as  the  response  time  to 
a  "DL/I  call"  in  an  IMS  data  base  management  system. 

We  consider  passage  time  simulation  of  this  network  under  the  following 
probabilistic  assumptions: 

(1)  for  1*0, 1,2,  and  3,  successive  service  times  form  a  sequence 
of  i.l.d.  random  variables; 

(11)  successive  £  service  times  form  a  sequence  of  i.l.d.  random 


variables; 
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(ill)  the  branching  at  successive  (respectively  a^)  service 

completions  is  governed  by  e  sequence  of  i.i.d.  Bernoulli  random 
variables  tp^  (respectively  tp2) » 

(iv)  the  sequences  in  (i) ,  (ii) ,  and  (iii)  are  mutually  independent; 

(v)  the  random  variables  aQ,  a^,  a2,  and  8  are  exponentially 
distributed; 

(vl)  P{ip2*l}>0,  so  that  jobs  eventually  receive  service. 

• 

The  circled  numbers  in  Figure  5.4  are  classes  of  jobs  in  the  indicated 
queues  and  are  the  key  to  representing  the  DL/I  component  model  as  a 
queueing  network  of  the  kind  described  in  Section  3.0.  This  representation 
is  shown  in  Figure  5.5.  Classes  2-7  are  identified  with  the  a  center  server 
(processor)  and  class  1  with  the  8  center  server  (I/O) .  At  the  a  center 
server,  classes  2-7  are  ordered  according  to  priority.  Class  7  service 
is  lnterruptable,  but  all  other  classes  of  service  are  not  interruptable. 

In  the  DL/I  component  model,  for  t2:0  we  define  the  vector  Z(t)  by 

Z(t)  -  (Q(t),C7(t),...,C2(t),S(t))  ,  (5.2.1) 

where 

Q(t)  -  number  of  jobs  waiting  or  in  service  at  8  center  at  time  t, 

S(t)  ■  class  of  job  being  served  by  a  center  server  at  time  t, 
and  for  231S7, 

C^(t)  ■  number  of  jobs  of  class  i  waiting  for  service  at  a  center 
at  time  t. 

We  let  N(t)  denote  the  position  at  time  t  of  the  marked  job  in  the  job 
stack  corresponding  to  the  order  of  the  nonzero  components  of  Z(t).  Then 


6. 


s»7 


Figure  5.5,  OL  / 1  component  model  with  job  classes 
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the  process  £*{X(t) :tkO},  where 

X(t)  -  (Z(t),N(t))  ,  (5.2.2) 

is  the  underlying  continuous  time  Markov  chain. 


The  Markov  chain  governing  change  of  job  classes  at  service  completions 
can  be  taken  to  have  state  space  (1,2 . 7}  and  transition  matrix  by 


Px  0  l-px  0 
0  l-p2  0  0 

px  0  l-px  0 
0  l-p2  0  0 

0  0  0  1 
0  0  0  0 


1 

0 

0 

0 

0 

0 

0 


0 

0 

0 

0 

0 

0 

1 


0 

0 

P2 

0 

p2 

0 

0 


where  p^  [respectively  p2]  denotes  the  probability  that  the  binary 
branching  variable  ^  [respectively  takes  the  value  1. 


The  sets  and  A2  defining  the  start  of  response  times  for  the  marked 
job  are 

A1  -  {(i,N-(i+l) ,0,0, 0,0, 0,7 ,N) :0£i<N> 

and 

A2  -  {(i,N-(i+l),0,0,0,0,0,6,N):Osi<N}  . 

Since  R  is  a  response  time,  Tj»Sj  for  jal,  and  B^*A^  and  B2~A2. 

Note  that  if  we  modify  the  DL/1  component  model  so  that  there  is  no 
preemption  of  service,  then  the  sets  A^  end  A2  are 
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Aj  -  {(i,N-(i+l),0,0,0,0,0,7,N)  rOSKN} 

and 

A2  -  { (i-m,N-(i+l) ,l,m-l,0,0,0,5,N-m) :l£i<N;  lSmSi} 
u  {(i,N-(i+l) ,0, 0,0, 0,0, 6, N) :0£i<N}  . 

This  provides  an  example  of  the  desirability  of  allowing  k<n  In  the  formal 
definitions  of  the  times  and  given  In  Section  4.2. 

5.3.  Numerical  Results 

The  results  displayed  here  were  obtained  using  the  particular  linear 
congruentlal  uniform  random  number  generator  described  by  Lewis.  Goodman 
and  Miller  (1969) .  Exponential  service  times  were  generated  by  logarithmic 
transformation  of  the  uniform  random  numbers;  Independent  streams  of 
exponential  random  numbers  (resulting  from  different  seeds  of  the  uniform 
random  number  generator)  were  used  to  generate  individual  exponential 
holding  time  sequences. 

The  numerical  values  of  point  estimates  and  confidence  intervals 
reported  here  for  r(f)»E{Y^(f ) }/E{a^}  are  the  classical  ratio  estimates, 
l.e.,  use  the  point  estimator  ?R(f)  of  Equation  (4.1.3)  and  the  100(l-2y)Z 

A 

confidence  interval  IR(f)  of  Equation  (4.1.4). 

Results  obtained  for  the  response  time  R  by  simulation  of  the  closed 
network  of  queues  of  Figure  5.1  are  in  Table  5.1.  The  initial  state  for 
the  Markov  chain  g  (and  return  state  identifying  cycles)  is  the  state 
(1,1).  The  results  in  Table  5.1  are  for  N«2  jobs,  with  u^-l,  y,®*0.5  and 
pa0.75.  Theoretical  values  for  the  long-run  expected  fractions  of  time 
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that  centers  1  and  2  ara  busy  (utilizations) ,  and  tha  expected  response 
time  can  ba  obtained  by  birth-death  process  methods.  These  values  ara  in 
parentheses.  In  all  cases,  the  90Z  confidence  intervals  obtained  surround 
the  theoretical  value. 

Results  of  simulation  for  characteristics  of  the  passage  time  P  are 
in  Tables  5.2  and  5.3.  For  N"2  jobs,  with  Uq"1»  Uj“0.5  and  <f>“0.75,  these 
tables  give  an  indication  of  the  effect  of  different  return  states.  The 
results  are  for  two  simulation  runs,  using  the  return  states  (1,1)  and 
(2,1),  respectively.  Thus,  for  example,  2367  transitions  in  the  Harkov 
chain  were  required  for  100  cycles  of  returns  to  the  state  (1,1)  but  only 
1183  transitions  for  100  cycles  of  returns  to  the  state  (2,1).  Since 
returns  to  the  state  (2,1)  occur  approximately  twice  as  frequently  as  to 
the  state  (1,1),  we  would  expect  that  only  half  as  many  cycles  (as  for 
return  state  (1,1))  for  comparable  accuracy;  this  Is  borne  out  by  the 
results  in  Tables  2  and  3.  In  all  cases,  the  90S  confidence  intervals 
obtained  surround  the  theoretical  value. 

Results  obtained  for  the  response  time  R  by  simulation  of  the  DL/I 
component  model  appear  in  Table  5.4.  The  initial  state  for  the  Markov 
chain  X  is  the  state  (1,N-2,0,0,0,0,0,6,N) .  Theoretical  values  for  the 
a  and  8  center  utilizations  along  with  the  expected  response  time  are 
obtainable  from  the  analysis  of  the  DL/I  component  model  given  by  Lavenberg 
and  Shedler  (1976).  Comparison  of  Table  5.4  with  Table  5.1  reveals  the 
effect  on  simulation  running  time  of  the  considerable  structural  complexity 
of  the  DL/I  component  model;  for  the  return  states  chosen,  there  are 
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i 

approximately  the  same  number  of  transitions  in  the  Markov  chain  for 
250  cycles  of  the  DL/I  component  model  as  for  400  cycles  of  the  simple 
model.  As  before,  the  90%  confidence  intervals  for  the  expected  response 
time  surround  the  theoretical  values. 


Slaulation  Results  for  Paaaaga  Tima  P  in  Cloaad  Network  of  Quauaa 
N-2,  p  -1 ,  W.-0.5,  p-0.75.  Return  State  la  (2,1). 


Slaulatlon  Results  for  Response  Tlae  R  In  DL/I  Coaponent  Model 
13,  E(a,)>l.S,  E{cu)-1,  ElBJ-50,  p  -0.1,  p,-0.2.  Return  State 
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6.0.  FINITE  CAPACITY  OPEN  NETWORKS  OF  QUEUES 

We  now  consider  open  networks  of  service  centers  with  jobs  srriving 
at  the  network,  traversing  the  network  and  receiving  service  along  the 
way,  and  finally  departing  from  the  network.  The  network  structure  we 
permit  is  essentially  the  same  as  that  described  in  Section  3.0,  except 
that  here  the  networks  are  open;  thus  arrivals  from  an  external  source 
and  departures  to  an  external  sink  occur.  Only  a  finite  number  of  jobs, 

N,  may  reside  in  the  network  at  a  given  time.  We  consider  two  formulations 
of  finite  capacity  open  networks.  In  "arrival  process  shutdown"  models, 
a  job  arriving  when  the  network  already  contains  N-l  customers  causes  the 
arrival  process  to  shut  down;  the  arrival  process  remains  shut  down  until 
the  first  subsequent  departure  from  the  network.  In  "jobs  turned  away" 
models,  the  arrival  process  never  shuts  down,  but  jobs  arriving  when  the 
network  already  contains  N  jobs  turn  away. 


6.1.  Markov  Arrival  Processes 

The  arrival  processes  we  allow  are  particular  stochastic  point 
processes  (series  of  events)  associated  with  time-homogeneous  Markov 
renewal  processes.  Let  J  be  a  finite  or  countable  set,  W*{Wn:n^0}  random 
variables  assuming  values  in  J,  and  U»{Un:nkO)  nonnegative  random  variables 
such  that  0-UqS;U1sU2£.  •  •  •  Recall  that  a  stochastic  process 
(W,U)-{(W  ,U  ) :nk0}  is  a  Markov  renewal  process  provided  that 


•  VrV'K . V  uo . 


■  p<VrJ-  WV'IV 


with  probability  one  for  all  nkO,  jeJ,  and  t£Q,  and  is  time-homogeneous 
provided  that  P{WQ+^«j ,  U^^-U^t  |Wq}  is  independent  of  n.  We  denote  by 
M»{M(t)  :t^0}  the  semi-Markov  process  generated  by  (W,JJ) ,  i.e., 

M(t)  -  W  .  if  U  s  t  <  U  .  .  (6.1.1) 

n  n  n+i 

By  definition,  a  Markov  arrival  process  is  a  stochastic  point  process  £ 
that  satisfies  the  condition 

..-U  sitlw  -i)  «  1  -  exp(-X.t)  (6.1.2) 

nri  n  n  l 

for  all  ieJ  and  t^O,  where  X^O.  Under  this  condition  the  semi-Markov, 
process  M  is  equivalent  to  a  continuous  time  Markov  chain.  Throughout, 
we  assume  that  M  is  irreducible  and  positive  recurrent,  and  that  the  Markov 
renewal  process  (W,U)  is  independent  of  the  service  times  and  Markov 
routing  of  jobs  in  the  network.  Note  that  in  "arrival  process  shutdown" 
models  we  can  think  of  the  arrival  process  as  operating  in  virtual  time, 
i.e.,  the  nth  customer  arrives  at  virtual  time  U  .  The  actual  time  of 

U 

the  nth  arrival,  however,  may  be  somewhat  later  due  to  the  finite  capacity 
constraint.  In  "Jobs  turned  away"  models,  the  nth  job  arrives  at  time 
U^,  but  may  be  turned  away. 

(6.1.3)  EXAMPLE.  Poisson  Process.  Let  J"{1}  so  that  W^-l  with  probability 
one  for  all  n20,  and  for  tiO  let 

p{VrVt}  ■ 1 "  exp(_Xt>  » 

where  X>0.  Clearly,  the  Markov  arrival  process  U«{U  :n£0)  is  a  Poisson 


process. 
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(6.1.4)  EXAMPLE.  Switching  Poisson  Process.  Let  J-{l,2, . . . ,k} , 
q_Lj»P{Wa^i»j  |Wu“l)  for  i,jcJ,  and  for  tkO  let 

p{VrVc|V1)  ■ 1  -  • 

where  X^>0  for  i«l,2,...,k.  Here  the  successive  times-between-events  ere 
exponentially  distributed  with  parameters  which  are  governed  by  the 
transition  matrix  (q^  :i,jcJ}.  This  Markov  arrival  process  Is  a 
semi-Markov  generated  point  process  with  exponential  holding  times. 

(6.1.5)  EXAMPLE.  Branching  Poisson  Process.  The  branching  Poisson  process 
is  a  model  for  clustered  arrivals.  Consider  a  particular  branching  Poisson 
process  constructed  as  follows.  A  Poisson  process  with  parameter  A,s 
generates  a  series  of  primary  events,  and  with  independent  probability 
rc(0,l],  a  primary  event  initiates  a  subsidiary  series  of  events.  Each 
subsidiary  process  consists  of  a  geometrically  distributed  number  of 
subsidiary  events,  and  the  times  between  these  subsidiary  events  are 
independent  and  exponentially  distributed  with  parameter  X^.  Finally, 

the  branching  Poisson  process  is  the  superposition  of  all  primary  and 

subsidiary  events.  To  represent  this  process  as  a  Markov  arrival  process, 

we  set  J*{0,1,2,. . .}  and  identify  W  with  the  number  of  subsidiary 

n 

processes  active  at  the  time  of  the  nth  event  of  the  branching  Poisson 
process.  Let  pc (0,1)  be  the  parameter  associated  with  the  geometric 
distribution  (mean  l/(l-p))  governing  the  number  of  subsidiary  events  in 

each  subsidiary  process;  that  is,  the  probability  of  k  subsidiary  events 

k-1 

in  a  given  subsidiary  process  is  p  q  (k2l) ,  where  q»l-p.  Then  for  ieJ, 


ft  \rrytw-*  y,  _ 


73 


!iA2q/(X1+i.X2)  ,  j-i-1 

IiX2p+(l-r)X1]/(X1+lX2),  j-i 
rX1/(X1+iX2)  ,  j-1+1  , 

and  for  t£0, 

p{VrVtfVi}  ■ 1  ■  «xpt-(ViX2)t3  • 

6.2.  Networks  of  Queues  and  Associated  Stochastic  Processes 

As  before,  we  permit  a  finite  number  of  single-server  service  centers, 
s,  and  a  finite  number  of  (mutually  exclusive)  job  classes,  c.  At  every 
time  epoch  each  job  Is  In  exactly  one  job  class,  but  jobs  may  change  class 
as  they  traverse  the  network.  Upon  completion  of  service  at  center  1,  a 
job  of  class  j  goes  to  center  k  and  changes  to  class  £  with  probability 
pij  kJL*  ,k£aYand  lSj.Xic)  and  £*{pj.  j^).  A  job  from  the  external 

source  arrives  at  center  k  as  a  job  of  class  £  with  probability  p^,  and 
a  job  of  class  £  completing  service  at  center  k  departs  to  the  external 
sink  with  probability  q^. 

At  each  service  center  jobs  queue  and  receive  service  according  to  a 
fixed  priority  scheme  among  classes;  the  priority  scheme  may  differ  from 
center  to  center.  Within  a  class  at  a  center,  jobs  receive  service 
according  to  a  fixed  discipline,  and,  some  centers  may  never  see  jobs  of 
certain  classes  as  determined  by  the  routing  matrix  P..  A  job  in  service 
at  a  center  may  or  may  not  be  preempted  if  another  job  of  higher  priority 
joins  the  queue  at  the  center.  All  service  times  in  the  network  are 
mutually  Independent,  and  at  a  particular  center  have  a  distribution 


i 
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consisting  of  exponential  stages  with  parameters  which  may  depend  on  the 
service  center,  class  of  Job  being  served,  and  "state”  of  the  entire 
network. 

Aa  in  Section  3.2,  we  let  S^(t) ,  i-l,2,...,s,  denote  the  class  of  job 

receiving  service  at  center  i  at  time  t,  with  S^(t)-0  if  at  time  t  there 

are  no  jobs  at  center  1.  The  classes  of  jobs  serviced  at  center  1  ordered 

by  decreasing  priority  are  J^iJ^Ci) . ^k(i)^’  el®ni®nt8  ot  th® 

set  {1,2 . c}.  Denote  by  C.^  (t)  ,C.^  (t) . (t)  the  number  of 

J1  J2  Jk(i) 

jobs  in  queue  at  time  t  of  the  various  classes  of  jobs  served  at  center  i. 

To  deal  with  passage  times,  we  again  use  the  idea  of  the  position  of 
a  distinguished  "marked"  job.  We  continue  to  think  of  the  (at  most  N) 
jobs  in  the  network  ordered  in  a  linear  stack  according  to  the  following 
scheme,  and  define  the  vector  Z(t)  by 

Z(t)  -  (C.(1)  (t),C.(1)  (t) . C.(1)(t),S  (t);...; 

Jk(l)  Jk(l)-1  J1 

C.°°  (t),C.(8)  (t),...,C.(8)(t),S  (t))  . 

Jk(s)  Jk(s)-1  J1 

The  linear  job  stack  corresponds  to  the  order  of  components  in  the 
vector  Z(t)  after  ignoring  any  zero  components.  Within  a  class  at  a 
particular  service  center,  jobs  waiting  appear  in  the  Job  stack  In  FCFS 
order,  the  latest  to  arrive  being  closest  to  the  top  of  the  stack.  We 
again  denote  by  N(t)  the  position  (from  the  top)  of  the  marked  job  in  this 
job  stack  at  time  t. 
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Recalling  that  M  is  the  semi-Markov  process  of  Equation  (6.1.1) 
associated  with  the  Merkov  arrival  process,  for  t£Q  we  specify  the 
state  of  the  network  by 

X(t)  -  (M(t),Z(t),N(t),U(t))  . 

Here  U(t)«(U, (t) , . . . ,U„(t))  and  U.  (t)  (lSjSN)  is  the  state  of  service 
1  N  j 

associated  with  the  jth  job-  in  the~job  stack;  (t)*Q  if  there 

are  less  than  j  jobs  in  the  system  at  time  t.  By  virtue  of  the  service 
time  distributional  assumptions,  the  Markovian  routing  structure,  and  the 
definition  of  a  Markov  arrival  process,  X"{X(t) : t£0}  is  a  continuous  time 
Markov  chain  with  a  (possibly  countable)  state  space  E. 

6.3.  Job  Marking 

The  principal  concern  here  remains  the  estimation  of  general 
characteristics  of  passage  times,  the  times  required  for  a  job  to  traverse 
a  specified  portion  of  the  open  network.  In  a  finite  capacity  open 
network,  a  passage  time  is  termed  a  "response  time"  if  it  is  the  total 
time  a  job  is  in  the  network.  To  estimate  passage  times,  we  track  an 
appropriate  sequence  of  typical  jobs,  based  on  the  idea  of  a  marked  job, 
and  measure  the  passage  times  for  a  sequence  of  marked  jobs;  these  are  to 
be  typical  jobs  in  the  sense  that  the  sequence  of  passage  times  for  the 
marked  jobs  should  converge  in  distribution  to  the  same  random  variables 
as  do  the  passage  times  for  all  jobs.  It  is  necessary  to  take  some  care 
to  ensure  that  this  is  the  case. 
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Our  lob  narking  schame  la  as  follows.  Let  i  (0Si<N)  be  the  number  of 
jobs  left  behind  by  the  marked  job  at  the  epoch  at  which  It  departs  from 
the  network  (i.e. ,  goes  to  the  external  sink).  Then  for  "arrival  process 
shutdown"  models,  the  (N-i)th  subsequent  arrival  is  the  next  marked  job. 
For  "jobs  turned  away"  models,  the  (N+l-i) th  subsequent  arrival  is  the 
next  marked  job. 

Note  that  this  marking  scheme  ensures  that  there  Is  at  most  one  marked 
job  in  the  network  at  a  time,  and  thus  there  is  no  need  for  further 
augmentation  of  the  numbers-in-queue ,  stages-of-service  state  vector  for 
the  measurement  of  passage  times.  (If  there  is  no  marked  job  in  the 
network  at  time  t,  we  define  the  position  N(t)  in  the  linear  job  stack  to 
be  zero.)  Note  also  that  more  than  one  passage  time  can  start  (and 
terminate)  for  a  particular  marked  job  before  it  departs  from  the  network. 

To  see  chat  the  sequence  of  passage  times  for  jobs  marked  by  this 
scheme  has  the  desired  property,  consider  "arrival  process  shutdown" 
models.  We  introduce  a  so-called  "phantom  server"  which  generates  the 
times  U*{UQ:niO}  according  to  the  Markov  renewal  process  (W,U)  of  the 
Markov  arrival  process.  Assuming  that  1  (0£1<N)  is  the  number  of  jobs 
left  behind  by  the  marked  job  at  the  instant  it  departs  from  the  network, 
we  route  the  marked  job  to  the  queue  at  the  phantom  server  where  N-l-1 
jobs  are  already  residing.  Upon  completion  of  service  to  a  job  by  tb ^ 
phantom  server,  the  job  returns  to  the  network  in  the  same  manner  as 
arriving  jobs,  i.e.,  with  probability  p^  the  job  goes  to  center  k  and 
becomes  class  l.  This  method  generates  arrivals  to  the  network  in  exactly 
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same  way  Che  original  arrival  process  does  with  Che  f inice  capaclcy 
conscrainC,  namely,  arrivals  occur  ac  Che  Clmes  UQ  provided  Che  number  in 
Che  necwork  is  less  chan  N.  Figure  6.1  illuscraces  Che  general  flow  of 
jobs.  Ic  is  now  clear  chac  in  effecC  we  have  a  closed  necwork  of  queues 
in  which  Che  marked  job  never  leaves  Che  necwork.  Furchermore,  Che 
scochascic  structure  of  Che  original  problem  remains.  The  chief  advancage 
of  chis  device  is  ChaC  we  can  use  Che  simulation  method  developed  in 
Secclon  3  for  closed  networks .  Note  Chat  there  is  no  need  for  a  "p<l" 
condition  Co  guarantee  stability  of  Che  open  network  since  ic  is 
effectively  a  closed  necwork  with  a  finite  number  of  jobs. 

For  "jobs  turned  away"  models,  we  introduce  a  phantom  server  in  a 
similar  way;  see  Figure  6.2.  A  job  completing  service  at  the  phantom 
server  with  j-1  other  jobs  ac  Che  server  returns  to  the  necwork  with 
probability  l-p(j),  where  p(j)«l  if  j»l  and  0  otherwise.  We  are  in  effect 
considering  a  closed  network  with  N+l  jobs.  It  is  straightforward  to 
check  that  the  simulation  method  of  Section  4.1  extends  to  closed  networks 
in  which,  as  here,  routing  probabilities  may  depend  on  the  number  of  jobs 
in  the  service  center. 

Having  reduced  the  problem  to  estimation  in  a  closed  network,  it  is 
necessary  to  modify  the  Markov  chain  X.  We  view  the  phantom  server  as 
service  center  s+1 ,  serving  only  one  class  of  jobs.  Let  Qs+^(t)  denote 
the  total  number  of  jobs  at  center  s+1  at  time  t.  Then  we  augment  the 
vector  Z(t)  with  the  component  Qg+^Ct),  and  define  X(t)  as  before  but  with 
this  augmented  Z(t) .  We  also  modify  the  Markov  routing  matrix  P»{p. .  . 

lj  jlCX 
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i  jobs 


Figure  6.2.  Flow  of  jobs  with  phantom  server  added. 
"Jobs  turned  away"  model 
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todescribe  the  dosed  network*  and  assume  that  the  resulting  matrix  is 
irreducible.  The  assumptions  made  on  the  network  imply  that  the  Markov 
chain  X  is  Irreducible  and  positive  recurrent,  and  it  is  possible  to 
proceed  with  the  estimation  of  passage  times  as  before. 

6.4.  Example  and  Numerical  Results 

We  consider  the  finite  capacity  open  network  of  queues  shown  in 
Figure  6.3.  Upon  completion  of  a  service  at  center  1  which  renders  a 
service,  in  accordance  with  a  binary- valued  variable  <(*,  the  job  joins  the 
tail  of  the  queue  in  center  1  (when  ^al)  or  joins  the  tail  of  the  queue 
in  center  2  which  renders  0  service  (when  .  Neither  center  1  nor 
center  2  service  is  subject  to  Interruption.  We  assume  a  FCFS  service 
discipline  in  each  queue.  The  limiting  response  time  of  interest  (denoted 
by  R)  is  the  time  measured  from  arrival  of  a  job  at  the  center  1  queue 
until  departure  of  the  job.  Also  of  possible  Interest  In  this  model  is 
the  limiting  passage  time  (denoted  by  P)  defined  as  the  time  measured  from 
arrival  of  a  job  into  the  center  1  queue  until  the  job  enters  the  center 
2  queue. 

For  each  of  the  arrival  processes  described  below,  we  consider 
simulation  of  this  network  of  queues  under  the  following  probabilistic 
assumptions: 

(1)  successive  ct  service  times  form  a  sequence  of  l.i.d.  random 
variables,  exponentially  distributed  with  rate  parameter  y^; 

(ii)  successive  6  service  times  form  a  sequence  of  l.i.d.  random 
variables  exponentially  distributed  with  rate  parameter  y^; 
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(i)  Service  center  1  renders  a  service 

(ii)  Service  center  2  renders  0  service 

(iii)  a  and  0  services  are  not  interruptabie 

(iv)  Routing  determined  by  binary  valued  random  variable 


Figure  6.3.  Finite  capacity  open  network  of  queues 
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(Hi)  is  a  Bernoulli  random  variable  and  value*  of  tj»  at  successive 

a  service  completions  form  a  sequence  of  i.i.d.  random  variables; 
(iv)  the  sequences  in  (i) ,  (ii)  and  (iii)  are  mutually  Independent. 

For  the  simulation,  the  generator  of  arrivals  is  either 

(i)  a  Poisson  process  (as  in  Example  (6.1.3))  having  rate  parameter  X; 
or 

(ii)  a  switching  Poisson  process  (as  in  Example  (6.1.4))  having 

parameters  A^,  A^,  p^  and  p2*  Here  X^  is  the  rate  parameter  of 
the  exponential  holding  time  in  state  1,  and  is  the  probability 
of  a  one-step  transition  from  state  i  to  state  i  in  the"embedded 
Markov  chain. 

In  this  network,  there  are  two  classes  of  jobs,  class  1  jobs  at  center 
1  and  class  2  jobs  at  center  2.  Since  each  center  sees  only  one  class  of 
job,  we  can  define 

Z(t)  -  (Q1(t),Q2(t))  , 

where 

Q^(t)  •  number  of  jobs  waiting  or  in  service  at  center  1  at  time  t  , 

and 

Q2(t)  •  number  of  jobs  waiting  or  in  service  at  center  2  at  time  t  . 

As  above,  let  N(t)  denote  the  position  of  the  marked  job  in  the 
job  stack  corresponding  to  the  order  of  the  nonzero  components  of  Z(t) , 
and  let  M(t)  denote  the  semi-Markov  process  associated  with  the  Markov 
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arrival  process-  Than,  sinca  the  service  times  are  exponentially 
distributed,  the  process  £"{X(t)-Tt20} ,  where 

X(t)  -  (M(t),Q1(t).Q2(t),N(t))  , 
is  the  underlying  continuous  time  Markov  chain. 

For  this  modal  the  state  space  E  of  X  is 

E  -  {(m.i.j.kJsOSi.J.kSN;  i+jSN;  m<sj}  , 
where  N  is  the  maximum  number  of  jobs  in  the  network. 

The  Markov  chain  governing  change  of  job  class  can  be  taken  to  have 
state  space  {1,2}  and  transition  matrix 


where  we  denote  by  p  the  probability  that  the  branching  variable  \j>  takes 
the  value  1.  The  probabilities  p^  governing  routing  of  jobs  from  the 
external  source  are  p^-1,  p12-0,  p2l"0,  p22*0,  and  tha  Probabilitie®  q^ 
of  departure  to  the  external  sink  are  q^-0,  q22*°»  q2l"0,  q22"1,  The 
3><3  matrix  governing  change  of  job  classes  for  the  associated  closed 
network  is 


' p  1-p  0 

0  0  1 

.1  0  0. 
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The  sets  and  defining  the  start  of  responae  tines  for  the  marked 
job  are 


A1  -  {(m,i,j,k)€E:k-N) 

and 

A2  -  {(n,i,j,l)6E:l>0>  . 

The  sets  and  defining  the  termination  of  response  times  for  the 
marked  job  are 


Bt.-  {(m,i,j,k)eE:k-i+j*  J>0} 

and 

B2  -  {(m,i,j,k)€E:k-i+j+l>  . 

Tables  6. 1-6. 3  contain  results  obtained  for  the  response  time  R  by 
simulation  of  the  network  of  Figure  6.3  for  N-4  jobs  with  Uq“1»  U^O.5 
and  p-0.75.  The  initial  state  for  the  Markov  chain  X  (and  return  state 
identifying  cycles)  is  the  state  (1,4, 0,1).  Simulation  results  for 
arrivals  generated  by  a  Poisson  process  of  rate  X"0.4  appear  in  Table  6.1. 
For  Poisson  arrivals  the  two  formulations  of  finite  capacity  networks  are 
equivalent.  Theoretical  values  for  the  long-run  expected  fractions  of 
time  that  centers  1  and  2  are  busy  (utilizations) ,  and  the  expected 
response  time,  obtainable  either  by  birth-death  or  semi-Markov  process 
methods,  appear  in  parentheses.  Note  that  with  the  exception  of  the  1000 
cycle  run,  the  90%  confidence  intervals  obtained  surround  the  theoretical 
value.  The  percentile  points  are,  respectively,  0.25,  0.5,  1,  1.5,  and 
2  times  the  theoretical  mean  response  time. 
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Table  6.2  gives  results  of  simulation  of  the  "arrival  process  shutdown" 
model  for  arrivals  generated  by  a  two-state  switching  Poisson  process  with 
X^«2,  \^m0.667,  Pj*0.978  and  P2*0.865.  With  these  parameter  values,  events 
in  the  stationary  switching  Poisson  process  occur  at  rate  0.4,  and  the 
stationary  tlmes-between-events  have  coefficient  of  variation  equal  to  2 
and  serial  correlation  coefficient  of  lag  1  (cf.  Cox  and  Lewis  (1966), 
p.  196)  equal  to  0.375.  In  all  cases,  the  90%  confidence  Intervals  for 
E{rJ  surround  the  theoretical  values.  Corresponding  results  for  the  "jobs 
turned  away"  model  are  in  Table  6.3.  Again,  the  90%  confidence  Intervals 
for  E{R}  surround  the  theoretical  value. 

An  overall  observation  from  Tables  6. 1-6. 3  is  that  the  lengths  of 
confidence  intervals  obtained  (for  expected  response  time  as  well  as 
percentiles  of  response  time)  from  the  three  simulations  are  roughly 
comparable. 
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10.0318  10.0246  10.0164  10.0111  10.0099  10.0084 


Simulation  Results  for  Response  Tima  R  In  Flnlta  Capacity  Open  Network  of  Queues. 
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*  I 

M  M 
<  -< 


P(R  s  23.836)  0.8760  0.6953  0.8979  0.8900  0.8898 

10.0218  10.0148  10.0121  10.0086  10.0078 


7.0.  MARKED  JOB  SIMULATION  VIA  HITTING  TIMES 


In  chis  section  we  develop  a  new  stochastic  setting  for  the  marlced 
job  method  based  on  "hitting  times"  of  an  underlying  continuous  time  Markov 
chain  to  fixed  sets  of  states.  This  formulation  is  the  basis  for  the 
subsequent  discussion.  We  first  make  some  remarks  about  the  results 
obtained  in  the  previous  sections. 

The  ealier  discussion  of  the  marked  job  method  applies  to  networks  of 
queues  having  single  server  service  centers.  The  generalization  to 
networks  having  multiple  server  service  centers,  however,  is 
straightforward.  To  handle  these  networks,  it  is  sufficient  to  incorporate 
into  the  job  stack  (and  associated  state  vector)  information  which  gives 
the  class  of  job  being  serviced  by  each  of  the  servers  at  a  multiple  server 
service  center. 

A  further  generalization  to  networks  with  stochastically  nonidentical 
jobs  is  also  possible.  The  marked  job  method  of  Section  4  applies  equally 
well  to  networks  in  which  jobs  are  of  a  finite  number  of  types.  In  such 
networks,  jobs  of  each  type  have  their  own  routing  structures  and  service 
requirements,  and  in  the  case  of  finite  capacity  open  networks.  Independent 
arrival  processes.  We  return  to  this  topic  in  Section  10. 

Although  the  finite  capacity  constraint  on  the  open  queueing  networks 
considered  in  Section  6  appears  natural  in  many  modelling  contexts,  it 
would  be  of  interest  to  extend  the  marked  job  method  to  open  networks  of 
infinite  capacity.  The  main  barrier  to  doing  so  is  the  definition  of  a 
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sequence  of  marked  jobs  whose  nonoverlapping  measured  passage  times 
converge  in  distribution  to  the  same  random  variable  as  does  the  entire 
sequence  of  passage  times  (enumerated,  say,  in  order  of  start  times). 

This  remains  an  open  problem. 

In  the  remainder  of  this  section,  we  develop  a  stochastic  setting  for 
passage  time  simulation  using  a  marked  job  which  is  simpler  than  that  of 
Section  4.2;  the  definition  of  passage  time,  however,  is  somewhat  more 
restrictive.  Here,  the  starts  of  passage  times  for  the  marked  job  are 
the  successive  entrances  (hitting  times)  of  a  particular  continuous  time 
Markov  chain  to  a  fixed  set  of  states,  and  the  terminations  of  such  passage 
times  are  the  successive  entrances  to  another  fixed  set  of  states.  This 
formulation  provides  the  basis  for  the  decomposition  method  for  passage 
time  simulation  in  Section  8,  the  analysis  of  statistical  efficiency  in 
Section  9,  and  the  discussion  in  Section  10  of  methods  for  networks  with 
multiple  job  types. 

7 . 1  Preliminaries 

As  in  Section  3,  we  consider  closed  networks  of  queues  with  a  finite 
number  of  jobs,  N.  In  each  network  there  are  a  finite  number  of  service 
centers,  s,  and  a  finite  number  of  job  classes,  c.  At  each  epoch  of  time 
each  job  is  in  exactly  one  job  class,  but  jobs  may  change  class  as  they 
traverse  the  network.  Upon  completion  of  service  at  center  i,  a  Job  of 
class  j  goes  to  center  k  and  changes  to  class  l  with  probability  p^ 

We  assume  that  P^tp^  j^lsi.kss,  lsj  ,JUc)  is  a  given  irreducible  Markov 
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AC  each  service  center  jobs  queue  and  receive  service  according  to  a 
fixed  priority  scheme  among  classes,  which  scheme  can  vary  from  center  to 
center.  Each  center  operates  as  a  single  server,  processing  jobs  of  a 
fixed  class  according  to  a  fixed  service  discipline.  All  service  times 
in  the  network  are  mutually  independent,  and  at  each  center  have  a 
distribution  with  a  Cox-phase  representation  with  parameters  which  may 
depend  on  the  service  center,  class  of  job  being  served,  and  the  "state" 
of  the  entire  system.  (Recall  that  we  exclude  zero  service  times  occurring 
with  positive  probability.)  A  job  in  service  may  or  may  not  be  preempted 
(according  to  a  fixed  procedure  for  each  center)  if  another  job  of  higher 
priority  joins  the  queue  at  the  center.  We  restrict  the  present  discussion 
to  networks  in  which  all  service  times  are  exponentially  distributed,  and 
deal  with  distributions  having  an  Cox-phase  representation  in  the  usual 
way  by  the  method  of  stages. 

We  review  the  notation  of  Section  3.2  used  to  characterize  the  state 

of  the  network  at  time  t.  For  1*1,2,.. .,s,  S^(t)  denotes  the  class  of  the 

job  receiving  service  at  center  1  at  time  t;  S^(t)-0  if  there  are  no  jobs 

at  center  i  at  time  t.  We  denote  by  j^U) , . . .  ,jk^  (i)  the  job  classes 

served  at  center  i,  ordered  by  decreasing  priority,  and  C^  (t) , . . .  ,C^  (t) 

J1  Jk(i) 

denote  the  number  of  jobs  in  queue  at  time  t  of  the  various  classes  served 
at  center  i. 

As  in  Section  3.2,  we  view  the  N  jobs  as  being  completely  ordered  in 
a  linear  job  stack,  and  define  the  vector  Z(t)  according  to 
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Z(t)  -  (C.(1)  (t),...,C.(1)(t),S.(t);...; 

Jk(l)  J1  1 

C.<8)  (t),...,C.(s)(t),S  (t))  .  (7.1.1) 

Jk(s)  31  8 

The  job  stack  corresponds  to  the  order  of  components  In  the  vector  Z(t) 
after  ignoring  any  zero  components.  Within  a  class  at  a  center,  Jobs 
waiting  appear  in  the  job  stack  in  the  order  of  their  arrival  at  the 
center,  the  latest  to  arrive  being  closest  to  the  top  of  the  stack. 

Letting  N(t)  denote  the  position  from  the  top  of  the  marked  job  in 
this  job  stack  at  time  t,  for  tiO  the  state  vector  of  the  network  is 

X(t)  -  (Z(t) ,N(t))  .  (7.1.2) 

As  before,  we  specify  the  passage  time  for  the  marked  job  by  four  subsets 

(A^,  A^,  B^,  and  B^)  of  E,  the  state  rpace  of  the  process  X«{X(t) :t£0). 

The  sets  A^  and  A^  [respectively  B^  and  B^]  determine  when  to  start 

[respectively  stop]  the  clock  measuring  a  particular  passage  time  for  the 

marked  job.  Denoting  the  jump  times  of  X  by  {t  :n£0},  for  k,n£l  we  require 

~  n 

that  the  sets  A^,  A2>  B^  and  B2  satisfy: 

U  X(To-l)eAl>  X(V4V  X(To-l+k)eAl  *nd  X<W4A2' 

then  X(t  . .  )«B.  and  X(t.  )eB.  for  some  0<msk; 
n-i*™  x  Ornn  £ 

and 

If  X<Vl>‘«l'  X<V‘*2'  XW!*1  *“>  X<Vk)eB2- 

then  X(TQ_^+a)eA^  and  X(tq+o)€A2  for  some  0*m<k. 

Recall  that  these  conditions  ensure  that  the  start  and  termination  times 
for  the  specified  passage  time  strictly  alternate.  Also  in  terms  of  these 
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jump  tines,  we  define  two  sequences  of  random  tines:  (S^ :J*0}  and 
{T^:j2l}.  The  start  [respectively  termination]  tine  of  the  Jth  pessage 
time  for  the  marked  job  is  denoted  by  [respectively  ] .  Assuming 
that  a  passage  time  for  the  marked  job  begins  at  t*0,  we  have 

SQ  -  0 

SJ  -  inf{Tn*Tj:X(Tn)«A2,  X(Tq_1)€A1},  j*l 

and 

Tj  -  inf{Tn>Sj_1:X(Tn)€B2,  X(tq-1)€B1},  jal  .  (7.1.3) 

The  jth  passage  time  for  the  marked  job  is  P  jfcl.  Note  that 

these  definitions  are  more  restrictive  than  those  of  Equation  (4.2.2)  in 
that  the  value  of  k  in  the  earlier  equation  is  set  equal  to  n-1. 

7.2.  Simulation  for  Passage  Times 

We  let  X  denote  the  state  of  the  continuous  time  Markov  chain  X  when 

U  'v 

the  (n+l)st  passage  time  of  the  marked  Job  begins:  X  «X(S  ),  n*0.  The 

n  n 

process  {X^n^Q}  is  a  discrete  time  Markov  chain  (with  state  space  A2) 
which  we  assume  to  be  irreducible  and  aperiodic.  By  Theorem  (4.2.5),  the 
process  { (Xq,  w  :n20)  is  a  regenerative  process  in  discrete  time,  and 
the  regenerative  property  guarantees  that  as 

<Xn*W  <X’P>  * 

The  random  variable  P  is  the  limiting  passage  time  for  the  marked  job. 
The  argument  in  Appendix  1  shows  that  the  sequence  of  passage  times  for 
any  other  job  also  converges  in  distribution  to  the  same  random  variable 
P.  Moreover,  the  sequence  of  passage  times  (irrespective  of  job  identity) 


94 


•numerated  according  to  start  tines  (or  termination  tinea)  also  converges 
to  P.  The  goal  of  the  simulation  is  the  estimation  of  characteristics  of 
the  limiting  passage  time  P.  Let  f  be  a  real-valued  measurable  function 
with  domain  R+,  and  define  the  quantity  r(f)  according  to 

r(f )  -  E{f(P)>  .  (7.2.1) 

We  now  depart  from  the  development  of  Section  4.2.  Let  L(t)  denote 
the  last  state  visited  by  the  Markov  chain  X  before  jumping  to  X(t),  end 
define  the  stochastic  process  V“{V(t) :t£0}  by 

V(t)  -  (L(t) ,X(t))  .  (7.2.2) 

In  this  development  of  the  marked  job  method ,  the  process  V  is  the 
fundamental  stochastic  process  of  the  passage  time  simulation. 

The  process  V  has  a  state  space  F  consisting  of  all  pairs  of  states 
(i,j),  i»JeE,  for  which  a  transition  in  X  from  state  1  to  state  j  can 
occur  with  positive  probability.  Since  J  is  an  irreducible,  positive 
recurrent  Markov  chain,  so  is  V.  The  entrance  times  of  V  to  a  state 
(l>j)eF  correspond  to  the  times  of  transition  in  £  from  state  1  to  state 
j.  We  define  subsets  S  and  T  of  F  according  to 

S  ■  {(k,m)«F:k«A1,  meAj) 

and 

T  -  { (k,m)eF:kcB^,  ,  (7.2.3) 

and  observe  that  the  entrances  of  V  to  S  [respectively  T]  correspond  to 
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the  starts  [respectively  terminations]  of  passage  times  for  the  marked 
job.  In  the  case  of  response  times,  S-T. 


The  next  step  is  to  select  a  fixed  element  of  S,  which  for  convenience 

we  designate  state  0,  and  assume  that  V(0)*Q.  We  let  {V^sniO}  denote  the 

embedded  jump  chain  associated  with  the  continuous  time  process  V.  The 

random  times  iy  :n£l}  denote  the  lengths  of  the  successive  0-cycles 
n 

(successive  returns  to  the  fixed  state  0)  for  {V  :n£0},  and  we  define  6“0 

u  U 

and  6  “Y,+. . ,+y  ,  mil.  Then  the  number  of  passage  times  for  the  marked 
ml  m 

job  in  the  first  0-cyde  of  V  is 


M 


1 


(Recall  that  for  a  set  A,  lA(x)-l  if  xeA  and  0  if  x^A.  Here  we  suppress 
the  argument  u.)  The  sum  of  the  values  of  the  function  f  for  the  passage 
times  for  the  marked  job  in  this  cycle  is  simply 


Yx(f) 


L f ' ‘v  • 

j-i 


We  denote  the  analogous  quantities  in  the  kth  0-cycle  by  and  Y^(f ) . 
Since  V  is  an  irreducible,  positive  recurrent  Markov  chain,  it  is  a 
regenerative  process,  and  the  pairs  of  random  variables  { (Y^(f) ,M^) :kil} 
are  Independent  and  identically  distributed.  Then,  provided  that 
P{PeD(f)}«0,  where  D(f)  is  the  set  of  discontinuities  of  the  function  f 
in  the  definition  of  r(f),  and  E{ |f (P) | }«*>,  it  follows  that 
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r(f)  -  E{f(P)}  -  E{Y1(f)}/E{M1>  . 

Therefore,  the  regenerative  method  applies  and  (from  a  fixed  rnanber  n  of 

0-cycles)  provides  the  strongly  consistent  point  estimate  fQ(f)-YQ(f)/Hn 

for  r(f).  Here  Y  (f)-(Y. (f)+. . .+Y  (f))/n  and  M  -(M.+...4M  )/n.  The 
n  i  n  ti  x  n 

associated  confidence  interval  is  based  on  the  central  limit  theorem 

nI/2{e„(f)-r(f)} 

n  ->  N(0,1)  » 

2 

where  o  is  the  variance  of  Y^(f)-r(f)M^.  It  is  easy  to  check  that  these 
point  and  interval  estimates  for  r(f)  (obtained  in  the  setting  of  the 
process  V)  are  the  same  as  those  obtained  from  Algorithm  (4.1.1). 

(7.2.4)  EXAMPLE.  Recall  the  model  of  Section  5.1  and  the  passage  time  P; 
see  Figure  7.1.  In  this  network  there  are  two  classes  of  jobs:  class  1 
jobs  at  center  1  and  class  2  jobs  at  center  2.  Since  each  center  sees  only 
one  job  class,  by  taking  into  account  the  fixed  number  of  jobs  in  the 
network,  for  t^O  we  can  define  Z(t)  to  be  the  number  of  jobs  waiting  or 
in  service  at  center  1  at  time  t.  Then  the  process  X-{(Z(t) ,N(t)) :tkO}, 
where  N(t)  is  the  position  of  the  marked  job  in  the  job  stack  at  time  t, 
has  state  space 

E  -  {(i,j) :0si£N,  lSJSN}  . 

For  the  passage  time  P,  the  sets  A^  and  A^  defining  the  starts  of 
passage  times  for  the  marked  job  are 
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^  -  {(i,N):OSi<H} 
and 

A2  -  {(i,l):0<isN>  .  (7.2.5) 

Similarly,  tha  secs  B^  and  defining  ehe  terminations  of  the  passage 
time  P  are 

Bx  -  {(i,i):0<iSN> 

and 

B2  -  {(i-l,i) :0<iSN}  .  (7.2.6) 

For  N-2  jobs.  Figure  7.2  shows  state  transitions  in  the  Markov  chain  X 
and  the  subsets  A^,  A2>  B^  and  B2  of  E.  The  process  V-{(L(t) ,X(t)) :tkO}, 
where  L(t)  is  the  last  state  visited  by  the  Markov  chain  X  before  jumping 
to  X(t) ,  has  state  space 

F  -  {(i,j,i+l,j+l):Osi<N,  lSj<N}  u  {(i,N,i+l,l) :0si<N}  u 
{(i,j,i-l,j):0<iSN,  UjSN)  u  {(i,i,i,l) :l<isN}  . 

The  subsets  of  F  defining  the  starts  and  terminations  of  passage  times 
for  the  marked  job  are 


S  -  {(i,N,i+l,l)  :0SKN} 


and 

T  -  {(i,i,i-l,i):0<isN)  . 


-’TV’?* 


gggBESgBgg 


(7.2.7) 
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8.0.  THE  DECOMPOSITION  METHOD 

In  this  section  we  concentrate  on  passage  times  through  a  subnetwork 
of  a  given  closed  network  of  queues,  i.e.,  passage  times  which  are  not 
complete  circuits  in  the  network.  For  such  passage  times,  we  develop  an 
estimation  method  in  which  observed  passage  times  for  all  the  jobs  enter 
into  the  construction  of  point  estimates  and  confidence  Intervals.  We 
consider  closed  networks  of  queues  as  in  Section  7.1.  The  formal 
specification  of  passage  times  is  as  in  Section  7.1,  but  we  make  the 
additional  assumption  with  respect  to  the  sets  S  and  T  of  Equation  (7.2.3) 
that 

S  a  I~  -  $  i 

this  effectively  rules  out  passage  times  which  are  complete  circuits,  i.e., 
response  times. 

The  basis  of  the  decomposition  method  for  estimation  of  passage  times 
through  a  subnetwork  is  simulation  of  the  network  in  random  blocks  defined 
by  the  terminations  of  certain  passage  times.  The  distinguished  passage 
times  are  those  that  (i)  terminate  when  no  other  passage  times  are  underway, 
and  (ii)  leave  a  fixed  configuration  of  the  job  stack  defined  by 
Equation  (7.1.1).  These  terminations  serve  to  decompose  the  sequence  of 
passage  times  for  all  of  the  jobs  into  Independent  and  identically 
distributed  blocks. 

We  denote  by  {P^:n2l}  the  sequence  of  passage  times  (irrespective  of 
n 

job  identity)  enumerated  in  order  of  passage  time  start.  As  before,  we 
let  f  be  a  real-valued  (measurable)  function  with  domain  R+,  and  the 
goal  of  the  simulation  is  the  estimation  of 
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r°(f)  -  E{f (P°) }  , 

where  P^>P^.  Note  that  by  the  results  of  Appendix  1,  P^-P ,  the  limiting 
n 

passage  time  for  (any)  marked  job. 

8.1.  Simulation  for  Passage  Times  Through  a  Subnetwork 

The  main  result  of  this  section  is  that  we  can  obtain  point  and 
interval  estimates  of  the  quantity  r^(f)  as  follows. 


(8.1.1)  ALGORITHM. 

1.  Select  a  configuration  z®  of  the  job  stack  at  which  a  passage 
time  terminates  and  there  are  no  other  passage  times  underway. 
Begin  the  simulation  with  this  configuration  of  the  job  stack  at 
time  0. 

2.  Carry  out  the  simulation  for  a  fixed  number  n  of  blocks  defined 
by  the  successive  terminations  of  passage  times  irrespective  of 
job  identity  which  leave  the  job  stack  in  the  (fixed) 
configuration  z®. 

3.  In  each  block,  measure  the  passage  times  for  all  the  jobs  and 
record  these  along  with  the  number  of  passage  times  observed  in 
the  block. 

4.  Denote  by  K®  the  number  of  passage  times  observed  in  the  mth 

m 

block  and  compute  Y^(f),  the  sum  of  the  quantities  f(p|?)  for  the 
m  J 

passage  times  Pj  in  the  mth  block,  e.g.. 


Y°(f) 
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5.  Take  as  a  point  estimate  (baaed  on  n  blocks)  of  r  (f)  the  quantity 

«“(£)  -  , 


where 


and 


-  a'1  £  *°<£> 


m*l 


K°  -  n'1  V  K°  . 
n  *  -t  m 


m*l 


6.  Take  as  a  100(1-2^)2  confidence  interval  (based  on  n  blocks)  for 
r^(f )  the  interval 


*n<f>  “ 


-0,-v  ,/-=0  l/2\  *0,-*  ,/-r0  l/2\ 

rn(f)'Zl-YSn/(V  /’  V^l-yVlV  / 


Here  sq  is  the  quantity 


sn  "  |911(a)-2f°(f)s12(n)+(r°(f)j2s22(n)|1/2  , 


where 


8ll(n)  -  (n-1) 


s22^  “  <n"1) 


-1  £«<*>■<)' 

S>“I 

m*l 


and 


’12 


(A)  -  to-l)'1^  (<<»<)(«) 


a-1/ 


The  quantity  (1-y) ,  where  $(•)  is  the  distribution  function  of  a 

standardized  normal  random  variable. 
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8.2.  The  Underlying  Stochastic  Structure 

We  begin  by  labelling  the  jobs  £roa  1  to  N,  and  for  1*1,2 , . . . ,N,  denote 
by  N^(t)  the  position  of  job  1  In  the  linear  job  stack  at  time  t.  Then, 
in  terms  of  the  vector  Z(t)  as  defined  In  Equation  (2.1.1),  we  set 

Xi(t)  -  (ZCc),!^)) 

and 

X°(t)  -  (Z(t)  y-Ct)  ,N2(t) . 1^(0)  .  (8.2.1) 

Each  of  the  processes  X**{X*(t) :tkO}  (1*1,2,. .. ,N)  and  g®»{X® (t) : tkO}  Is 
an  irreducible,  positive  recurrent  continuous  time  Markov  chain.  We  denote 
the  state  space  of  the  process  X1  [respectively  X®]  by  E1  [respectively  E^] . 

Next  we  let  L*(t)  [respectively  L^(t) ]  denote  the  last  state  visited 
by  the  Markov  chain  X1  [respectively  X°]  before  jumping  to  Xi(t) 

[ respectively:  X^ (t>]  J:  and  for  1*1,2, . .  _,N  and  t£0  define 

Vi(t)  -  (Li(t),Xi(t)) 

and 

V°(t)  -  (L°(t) ,X°(t))  .  (8.2.2) 

The  process  V^*(v®(t) :t*0}  is  the  fundamental  stochastic  process  of  the 
simulation  for  passage  times  through  a  subnetwork.  Note  that  incorporation 
of  the  component  L^(t)  into  the  definition  of  V®  is  necessary  for  detection 
of  the  starts  and  terminations  of  passage  times.  Since  each  of  the 
processes  X*  and  X®  is  an  irreducible,  positive  recurrent  continuous  time 
Markov  chain,  so  is  each  of  the  processes  V^’»{Vi(t)  :t*0}  and  V®.  We  denote 
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the  state  spaces  of  V*  and  by  F*  and  F^,  respectively.  As  In 
Equation  (7.2.3),  we  define  subsets  S*-  and  T*  of  F^  according  to 

S*  ■  { (k,m)«F*:k«A^,  m€A2} 

and 

T1  -  {(k,m)€Fi:kcB1,  me^}  .  (8.2.3) 

Thus,  the  entrances  of  V*  to  S*  [respectively  T*]  correspond  to  the  starts 
[respectively  terminations]  of  passage  times  for  job  1.  Ue  also  define 
two  subsets  S®  and  of  F®  according  to 


and 


{ (z,n^, . . . ,n^,z' ,n^, . . . ,n^)cF®:  for  some  k, 
(z,nk)eA1  and  (z',n£)eA2} 


{(z.n^, ..  .n^.z*  ,n^,. . .  ,n^)eF°:  for  some  k, 
(z.n^eBj^  and  (z',n£)eB2}  . 


(8.2.4) 


The  entrances  of  to  the  set  correspond  to  the  starts  of  passage 
times  (Irrespective  of  job  identity)  and  the  entrances  of  V®  to  the  set 
T®  correspond  to  the  terminations.  Thus,  from  a  simulation  of  the  process 
V  ,  it  Is  possible  to  measure  the  passage  times  for  each  of  the  jobs. 


Now  consider  {P^snkl}  the  sequence  of  passage  times  (irrespective  of 
n 

job  identity) ,  enumerated  in  order  of  passage  time  start.  Formal 
definition  of  the  sequence  {P®}  is  in  terms  of  sequences  of  starts  and 
terminations  of  passage  times  for  each  of  the  jobs;  the  definitions  of 
the  latter  involve  entrances  of  ^  to  the  sets  S*  and  T1  of  Equation  (8.2.3) 
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and  are  analogous  to  Equation  <7.1.3).  Recall  that  the  goal  o£  the 
simulation  Is 

r°(f)  -  E{f(P°)>  ,  (8.2.5) 

where  P^— »P^,  and  £  Is  a  real-valued  measurable  function  with  domain  R. . 

XI  T 

We  carry  out  the  simulation  In  random  blocks  o£  the  process  V®  defined 
by  the  successive  entrances  of  V®  to  a  fixed  set  of  states  0®.  Entrances 
of  V°  to  the  set  U°  correspond  to  the  terminations  of  passage  times 
(Irrespective  of  job  identity)  which  occur  when  no  other  passage  times 
are  underway,  and  which  leave  a  fixed  configuration  of  the  job  stack. 

Formally,  let  D  be  the  state  space  of  the  process  Z"{Z(t):tiO}  defined 
by  Equation  (3.2.1),  and  denote  by  C  the  set  of  (center,,  class)  pairs  In  the 
network.  We  define  a  function  h  taking  values  In  C  and  having  domain 
D*{l,2, . . . ,N}  as  follows.  For  zeD  and  ne{l,2 , . . . ,N> ,  the  value  of  h(z,n) 

Is  (l,j)  when  the  job  in  position  n  of  (the  job  stack)  z  is  of  class  j  at 
center  1.  Now  consider  the  embedded  jump  chain  {V^k^O}  associated  with 
the  continuous  time  Markov  chain  V  of  Equation  (7.2.2).  For  states  v' ,v”«F 
the  state  space  of  {V^:k20},  we  write  v'o^v"  when  v"  is  accessible  from  v' , 
i.e.,  when  for  some  nkl,  the  probability  starting  from  v'  of  entering  v" 
on  the  nth  step  is  positive.  Similarly,  for  any  subset  L  of  F  we  write 
v's^rv"  when  v"  is  accessible  from  v’  under  the  taboo  L;  this  means  (cf.  Chung 
(1967),  pp.45,  48)  that  for  some  nil,  there  is  a  positive  probability, 
starting  from  state  v',  of  entering  state  v"  on  the  nth  step  under  the 
restriction  that  none  of  the  states  in  L  is  entered  in  between  (exclusive 
of  both  ends) . 


*  '  .'A- 
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Next,  we  define  a  subset  H  of  C  according  to 

H  -  {(i,j)£C:  for  some  (z,n,z' ,a*)€T-S,  h(z’,n')  ■  (i,j)}  u 
{(i,j)eC:  for  some  (z,n,z' ,n')eF-(TuS) ,  v'cT  and  v"eS, 
v'/v*(z, n,z',n’),  (z,n,z'  ,n’)  &v"  and  h(z'  ,n’)  -  (l,j)} 

Thus,  the  set  where  a  (center,  class)  pair  is  in  [respectively 

H,]  if  it  is  possible  for  the  marked  job  to  be  of  class  j  at  center  1  when 
the  passage  time  specified  by  the  sets  ,  B^,  and  terminates 

[respectively  is  not  underway].  Note  that  the  set  H  is  nonempty  since  by 
assumption  SnT«$  and  thus  is  nonempty. 

Now  define  a  subset  of  D»  the  state  space  of  {Z(t):t£0},  according  to 

D°  -  (zeD:h(z,n)eH  for  n-l,2,,..,N  and  for  some  n,  Mz.nJcH^  . 

Elements  of  the  set  D°  correspond  to  configurations  of  the  job  stack  upon 
termination  of  a  passage  time  with  no  other  passage  times  underway.  The 
set  is  nonempty  since  H  is  nonempty.  Therefore,  we  can  select  an 
element  z^  of  D®,  and  in  terms  of  this  fixed  z®,  finally  define  the  set 
according  to 

U°  -  { (*,nx . nN,z,,n^,...,n^)6T°:z’-z0}  ,  (8.2.6) 

where  is  given  by  Equation  (8.2.4). 

For  convenience,  we  assume  that  V^(0)«U^.  The  random  times  {Y®:m£l} 

m 

denote  the  lengths  of  the  successive  blocks  (returns  to  the  set  U°)  for 
{V^:n20}  the  embedded  jump  chain  associated  with  V°,  and  we  define  6^-0 
and  6^*Y?+. . m2l. 

m  i  m 
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The  number  of  passage  times  in  tha  first  block  of  tha  process  V® 
is 

and  we  denote  the  analogous  quantity  in  the  mth  block  of  by  K^.  Mote 
that  within  each  block  of  V®  defined  by  the  entrances  to  the  set  ,  at 
least  one  passage  time  starts  and  terminates. 


Next,  we  let  Y°(f)  be  the  sum  of  the  quantities  f (P?)  over  the  passage 
m  J 

times  Pj  in  the  mth  block  of  V^,  for  example, 


Proposition  (8.2.7)  contains  the  key  observation. 


(8.2.7)  PROPOSITION.  The  sequence  of  pairs  of  random  variables 

{(Y°(f),K°) :m2l)  are  Independent  and  identically  distributed, 
m  m 

The  argument  used  in  Appendix  1  shows  that  P^-oP^  as  and  that 

this  random  variable  P°  is  the  same  random  variable  as  the  limiting  passage 
time  P  of  (any)  marked  job.  For  the  function  f  appearing  in  the  definition 
of  r°(f),  let  D(f)  denote  the  set  of  discontinuities  of  f.  Assuming  that 
P{P°cD(f)}«0  and  using  Lemma  (2.19),  it  follows  that  as  tr*»,  f (P®)“Of (P°) . 

U 

Finally,  standard  arguments  (cf.  Appendix  2)  yield  a  ratio  formula  for 
r°(f) . 
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(8.2.8)  PROPOSITION.  Assume  tbet  E{ | f (P°) |}<».  Then 

r°(f)  -  E{f (P°) }  -  E{Yj(f)}/E{Kj}  . 

With  Che  ratio  formula  of  Proposition  (8.2.8)  and  the  result 

(Proposition  (8.2.7))  that  the  pairs  of  random  variables  {(Y®(f) ,K^) :m£l) 

are  i.i.d.,  the  regenerative  method  applies;  from  a  fixed  number  of  blocks 

we  obtain  the  point  estimate  r°(f)-Y°(f)/K°  and  an  associated  confidence 

n  n  n 

interval  for  r®(f). 

(8.2.9)  EXAMPLE.  Consider  the  model  of  Section  5.1  and  Example  7.2.4. 

First  observe  that  for  the  passage  time  P,  the  sets  S  and  T  (given  by 
Equation  (7.2.7))  are  disjoint.  The  process  X°«{ (Z(t)  .^(t)  , . . .  .N^t))  :t*0} 
has  state  space  E°,  where 

E°  -  {(i,^,...,!^)  :0si£N;  1*^ . r^SN;  n^^  for  i^j}  . 

The  underlying  continuous  time  process  V®-{(L®(t) ,X^(t)) :tiO),  where  L^(t) 
is  the  last  state  visited  by  the  Markov  chain  X°  before  jumping  to  X°(t), 
has  state  space  F° .  The  subsets  of  F°  defining  starts  and  terminations 
of  passage  times  are 

o  * 

S  -  {(i,n1,...,nN,i+l,n^,...,n^l):  for  exactly  one  j  ,  n  *-N 

and  n\-l;  lSn^<N  and  nj-n^+l,  ji*j*;  n^n^  for  i#j} 

and 

T°  -  {(i,ni . nN,i-l,nJ>...,n^):0<iSN;  lsnj-njsN, 

for  lSjSN;  for  lt*j }  . 


•  i.  . 
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For  this  network  the  process  £«{Z(t) stkO)  hes  stete  spece  D«{0,1, . . . ,N>. 

The  set  C  of  (center,  class)  pairs  Is  {(1,1), (2 ,2)},  the  set  ^*4  end 
H<-Hj-{(2,2)}.  The  set  D°-{0},  and  since  D®  Is  a  singleton,  we  must  select 
2®»0.  For  z^*0,  the  set  U®  defining  the  blocks  of  the  process  V®  Is 

U°  -  {(l,n1,...,nN,0,n1,...,nN):lsn1 . n^sN;  n^^ ;  for  i*j}  . 

Now  consider  the  case  of  N«2  jobs.  The  state  space  F®  of  the  continuous 
time  Markov  chain  has  ten  elements,  and  the  set 

U°  -  {(1,1, 2, 0,1, 2),  (1,2, 1,0, 2,1)}  ; 

see  Figure  8.1.  A  portion  of  a  sample  path  for  appears  in  Figure  8.2, 
and  Figure  8.3  shows  the  corresponding  decomposition  of  the  sequence  {P^} . 

For  passage  times  through  a  subnetwork,  the  decomposition  method 
provides  an  alternative  to  marked  job  simulation.  Since  observed  pessage 
times  for  all  of  the  jobs  enter  in  the  construction  of  these  point  and 
interval  estimates,  we  would  expect  this  method  to  have  greater  statlsticel 
efficiency  than  the  marked  job  method.  In  this  connection,  the  celculatlon 
of  theoretical  values  for  variance  constants  entering  into  central  limit 
theorems  used  to  obtain  confidence  intervals  from  passage  time  simulations 
la  of  interest.  These  calculations  are  the  subject  of  Section  9. 


•  - '  ’  -^IKCWWI H  ■*&'*£  • 


S|_ ,  ^Start  time  for  jth  passage  time  for  job  i,  i-1 ,2 
T| -Termination  time  for  jth  passage  time  for  job  i,  i-1 ,2 
X  Entrances  of  V°  to  U° 


Figure  8.2.  Sample  path  for  process  V° 


Figure  8.3.  Decomposition  of  sequence  of  passage  times 


3 


9.0.  EFFICIENCY  OF  SIMULATION 

We  consider  in  this  section  the  calculation  of  theoretical  values  for 
variance  constants  entering  into  the  central  limit  theorems  used  in 
previous  sections  to  obtain  confidence  Intervals  for  passage  time 
characteristics.  Using  results  of  Hordijk,  Iglehart  and  Schassberger 
(1976)  for  the  calculation  of  moments  in  discrete  time  and  continuous  time 
Markov  chains,  we  compute  theoretical  values  for  mean  passage  times.  We 
do  this  first  for  the  marked  job  method  (in  the  stochastic  setting  of 
Section  7),  and  then  for  the  decomposition  method  of  Section  8.  For  passage 
times  where  both  estimation  methods  apply,  these  results  provide  a  firm 
basis  for  a  comparison  of  statistical  efficiency.  The  calculations  also 
make  it  possible  to  assess  the  efficacy  of  the  marked  job  method  for 
simulation  of  response  times. 

9.1.  Theoretical  Values  for  Finite  State  Markov  Chains 

Following  Hordijk,  Iglehart,  and  Schassberger  (1976),  we  first  consider 
discrete  time  Markov  chains,  and  let  {X^:k20}  be  an  Irreducible  Markov 
chain  with  finite  state  space  E*{0 ,1, . . . ,N}  and  one-step  transition  matrix 

l  -  {pi:j:i,jeE}  . 

For  this  chain,  we  let  denote  the  n-step  transition  probability  from 

state  1  to  state  j,  and  recall  that  for  n21, 

In  -  (p”j :i,j€E)  . 

Throughout  this  section  we  use  the  following  notation.  For  a  fixed 
state  ieE,  P^{ * }  denotes  the  conditional  probability  associated  with 


.HI.  1N1UH.  L»U1 ■ 
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starting  the  chain  In  state  1,  and  E^{  * }  denotes  the  corresponding 
conditional  expectation.  For  jeE,  the  state  space  o£  {X^:k20},  and  nkl 
we  let  $Q(j)  denote  the  nth  entrance  time  of  {X^:k20}  to  state  j,  e.g., 

0x(j)  -  mlntkiltX^-j} 

and  let  a^(j)«8j.(j)  and  a  (J)«8nCj)-B  ^(J) ,  n>l.  This  notation  is 

consistent  with  that  introduced  in  Section  2  for  regenerative  processes. 
Note  that  {$Q(j):nkl}  is  a  (possibly  delayed)  renewal  process  since  a 
finite  state,  irreducible  Markov  chain  is  nacessarily  positive  recurrent 
and  therefore  returns  to  every  state  jeE  infinitely  often  with  probability 
one.  If  Xq*J ,  the  process  {8Q(j):nkl}  is  an  ordinary  renewal  process. 

We  consider  vectors  such  as  (v(0) ,v(l) , . . . ,v(N))  to  be  column  vectors. 
Real-valued  functions,  such  as  f  and  g,  having  domain  E  are  viewed  in  this 
way  and  denoted  by  f_  and  &.  In  this  context  the  symbol  E{*}  denotes  the 
vector 


(e0{*},e1{-},...,en(-})  . 

In  addition  (for  vectors  u  and  v)  the  symbol  u°v  denotes  the  vector 

(u(O)v(O) ,u(l)v(l) .... ,u(N)v(N))  . 

For  a  matrix  A»(ag,a^, . . . ,aa> ,  we  let 

u°A  •  A®u  ■  (u°a_ ,u°a, , . . . ,u°a  )  . 

u  i  m 

Finally,  for  a  matrix  B«(bn ,b. , . . . ,b  ),  we  let 

“  u  x  m 

A®B  -  (a0ob0,a1»b1,...,aa«ba)  . 


i 
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We  first  consider  computation  of  the  values  of  E^{Y^(f)}  and 
EffYiCf )Yi(g) }  for  general  real-valued  functions  f  and  g  having  domain  E, 
and  ieE.  Accomplishing  this  (and  similar  computations  for  continuous  time 
Markov  chains)  makes  it  possible  to  obtain  theoretical  values  for  mean 
passage  times;  in  addition  we  can  obtain  values  for  the  variance  constants 
which  enter  into  central  limit  theorems  used  to  obtain  confidence  intervals 
in  simulations  by  the  marked  job  and  decomposition  methods. 


For  the  discrete  time  Markov  chain  (X^rkkO),  we  consider  here  only 

cycles  of  the  regenerative  process  formed  by  the  successive  entrances  to 

state  0,  and  henceforth  suppress  the  0  in  the  notation  B  (0) ,  a  (0),  etc. 

n  n 

Note  that  this  is  no  real  restriction,  and  that  equally  well  we  could 
choose  any  other  state  ieE.  For  i.jeE  and  n-0,1 .  let 


0Pij  *  Pi{cti>n»  Vj}  * 

and  set 

“  toPij:i,jeE}  * 

We  obtain 

m  p 

from  P  by  setting  the  0-column  of  P  equal  to  0.  It  is  easy  to  see  that 
qP11  is  the  matrix  product  of  n  copies  of  qP,  and  that  for  all  nkl, 


n 

0pi0 


0  . 


•>  1  * 
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Theorem  (9.1.1)  Is  due  to  Hordijk,  Iglehert,  end  Schessberger  (1976). 
For  any  real-valued  function  f  with  doaeln  E,  we  define 

ei_1 

Y  (f)  -  £  fOL)  • 

A  k-0 


(9.1.1)  THEOREM.  For  an  Irreducible,  finite  state  discrete  time  Markov 
chain  with  transition  matrix  P_, 

E{Y1(f)}  -  (I-qP)"1!  (9.1.2) 

and 

E{Y1(f)Y1(g)}  -  (!-(£) _1lL  ,  (9.1.3) 

where  h-f^EtY^  (g)  H£«E{Yj  (f  ) 


How  we  consider  continuous  time  Markov  chains  and  let  X-{X(t) :t£0}  be 
a  Markov  chain  with  finite  state  space  E«{0,1, . . . ,N} ,  transition 
matrix  F>(t)«{pij(t) :l,jeE}  and  matrix  of  infinitesimal  transition  parameters 
(^{q^ :i,jcE} .  Recall  that  in  a  continuous  time  Markov  chain  the 
matrix  <JvP' (0)  is  the  given  data.  In  general,  £(t)  is  difficult  to 
calculate  and  is  rarely  available  explicitly.  The  exponentially 
distributed  holding  time  in  any  state  icE  has  rate  parameter  q^-q^.  For 
all  leE,  we  assume  that  0<q^<°°,  i.e.,  that  all  states  are  stable  and 
nonabsorbing,  and  in  addition  chat 


N 


E 

J-0 


qiJ 


0  . 
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This  last  assumption  guarantees  that,  starting  from  any  state  leE,  the 
Markov  chain  X  makes  a  transition  to  a  next  state  jeE.  We  now  form  the 
jump  matrix  R*{r^j}  of  the  chain,  defining  the  elements  r^  according  to 

qij/qi  '  l*1 
0  ,  j-i 

We  assume  that  the  jump  matrix  R  Is  Irreducible  (and  therefore  positive 
recurrent) .  This  is  equivalent  to  the  continuous  time  Markov  chain  X  being 
irreducible.  As  in  the  case  of  a  discrete  time  Markov  chain,  we  let  P^{*} 
and  E^*}  denote  the  conditional  probability  and  conditional  expectation 
associated  with  starting  in  state  ieE.  For  jeE  and  n£l,  we  let  BQ(j)  denote 
the  nth  entrance  time  of  X  to  state  j,  i.e., 

^(j)  -  inf{s>0:X(s'-)^J ,  X(s)-j}  . 

We  now  consider  the  computation  of  E^{Y^(f)}  and  ^{Y^(f)Y^(g)  }  for 
real  valued  functions  f  and  g  with  domain  E.  As  in  the  case  of  discrete 
time  Markov  chains,  we  restrict  attention  to  regenerative  cycles  formed 
by  the  successive  entrances  to  state  0,  and  suppress  the  0  in  our  notation. 
For  t20,  we  let 

qP^U)  -  P1(a1>t,  X(t)-j)  , 

-  {^^(t)  :i,jeE}  , 

and,  for  n£0,  construct  the  matrix  ^Rn  from  R  in  the  same  manner  as  we 
constructed  gPn  from  P  in  the  discrete  time  case. 
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For  *  real-valued  function  f  defined  on  Et  we  define  Y^(f  )  according 
to 


Yx(f)  -  1  f (X(t))dt  , 

and  let  be  the  column  vector 

-1  ,  -1  -1  -1* 

1  “  W0  •*•••%  )  • 

Theorem  (9.1.4)  is  due  to  Hordijk,  Xglehart,  and  Schassberger  (1976). 


(9.1.4)  THEOREM.  For  an  irreducible,  finite  state  continuous  time  Markov 
chain  with  jump  matrix  R  and  vector  q  of  rate  parameters  for  holding  times. 


and 


E{Y1(f)>  -  e|J  o2(t>£dtJ  “  (I-qR)'1^"1) 
E{Y1(f)Y1(g)}  -  e|J^  ^P(t)h  dtj  -  (I-qR)"1^"1) 


(9.1.5) 


(9.1.6) 


where  h"f^*E{Y^(g)}+£,,E{Y^(f)}. 


He  now  show  how  to  use  the  results  of  Theorems  (9.1.1)  and  (9.1.4)  to 
assess  the  statistical  efficiency  of  simulation  by  the  marked  job  method 
for  mean  passage  times. 


9.2.  Variance  Constants  for  the  Marked  Job  Method 

We  consider  closed  networks  of  queues  and  passage  times  as  in 
Section  7.1.  For  t^O,  the  state  vector  of  the  network  is 


_ . _ _ r 1  ■-"» '  rn  ~ 
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X(t)  -  (Z(t).N(t))  , 


(9.2.1) 


where  Z(t)  (of  Equation  (7.1.1))  corresponds  to  the  linear  job  stack,  and 
N(t)  Is  the  position  In  the  job  stack  of  the  marked  job  at  time  t.  Recall 
that  the  process  X-{X(t) : t^O}  Is  an  Irreducible,  positive  recurrent  Markov 
chain  with  state  space  E.  As  before,  we  denote  by  L(t)  the  last  state 
visited  by  the  Markov  chain  X  before  jumping  to  X(t) ,  and  the  process 
V-{V(t):tiO}  defined  by 


V(t)  -  (L(t) ,X(t))  (9.2.2) 

Is  the  fundamental  stochastic  process  of  the  passage  time  simulation. 

Recall  that  the  process  V  has  a  state  space  F  consisting  of  all  pairs 
of  states  (l,j),  i,jc£  for  which  a  transition  in  X  from  state  1  to  state  j 
can  occur  with  positive  probability.  Since  X  Is  an  Irreducible,  positive 
recurrent  Markov  chain,  so  is  V,  and  the  entrances  of  V  to  the  fixed 
subset  S  [respectively  T]  (defined  by  Equation  (7.2.3))  of  the  state  space  F 
correspond  to  the  starts  [respectively  terminations]  of  passage  times  for 
the  marked  job. 

As  In  Section  7.2,  we  select  a  (fixed)  state  of  S,  designated  state  0, 
and  assume  that  V(0)*0.  To  estimate  the  quantity  r(f)  of  Equation  (7.2.1), 
the  marked  job  method  prescribes  that  we  carry  out  the  simulation  of  V  in 
0-cycles  defined  by  the  successive  returns  to  state  0;  within  each  cycle 
we  record  the  number  of  passage  times  of  the  marked  job  and  measure  each 
of  these  passage  times. 
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The  key  results  of  Section  7.2  lending  to  point  estimates  end 
confidence  Intervals  for  the  quantity  r(f)  are  that  the  pairs  of  random 
variables 


{(Yk(f),Mlt):k2l}  (9.2.3) 

are  Independent  and  Identically  distributed,  and  that 

r(f)  -  EQ{Y1(f) }/Eq{M1}  .  (9.2.4) 

Recall  that  is  the  number  of  passage  times  for  the  marked  job  In  the 
kth  0-cycle  and  Y^(f )  Is  the  sum  of  the  values  of  the  function  f  for  the 
passage  times  of  the  marked  job  in  this  cycle. 

Given  Equations  (9.2.3)  and  (9.2.4),  the  regenerative  method  provides 
(from  a  fixed  number  n  of  0-cycles)  the  point  estimate 

r(f)  -  Y  (f)/ff  . 
n  n  n 

The  associated  confidence  interval  for  r(f)  follows  from  the  central  limit 
theorem 

n1/2{?  (f)-r(f)} 

- S - -  N(0,1)  ,  (9.2.5) 

O(f)/E0{MX) 

where 

a2(f)  -  var{Y1(f)-r(f)M1)  .  (9.2.6) 

For  calculation  of  theoretical  values,  we  restrict  attention  to  the 
mean  passage  time;  thus,  the  function  f  in  the  definition  of  r(f)  is  the 
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Identity  function.  Using  the  results  of  Section  9.1,  we  show  how  to  compute 

the  value  of  the  mean  passage  time  r  and  the  corresponding  variance 
2 

constant  0  appearing  In  the  central  limit  theorem  of  Equation  (9.2.S). 

These  computations  rest  on  the  definition  of  two  particular  functions 

*  *  .  . 
(denoted  f  and  g  )  having  domain  F  and  taking  values  in  the  set  10,1/. 

* 

We  define  the  function  f  to  be  the  indicator  function,  lg,  of  the 
set  S  which  defines  the  starts  of  passage  times  for  the  marked  job;  i.e., 
for  (z,n,z' ,n')«F, 

f*(z,n,z' ,n')  -  lg(z,n,z' ,n)  .  (9.2.7) 

Proposition  (9.2.8)  follows  directly  from  Theorem  (9.1.1). 

(9.2.8)  PROPOSITION.  Let  f *  be  the  function  defined  by  Equation  (9.2.7), 
and  R  the  transition  matrix  of  the  discrete  time  Markov  chain  {V^:k20}. 

Then 

E{Y1(f*)}  -  E 
and 

E{(Yx(f*))2}  -  (I-qRTV  , 

where  6^  is  the  time  of  the  first  return  to  the  state  0  in  {V^:k£0}  and 

*  *  *  ,  *  * 
h  «2f  ®E{Y1(f  )}-f  of_  . 

We  use  Proposition  (9.2.8)  and  the  definition  of  to  obtain  the 

2 

quantities  EgtM^}  and  Eq{M*}  according  to 


fir1 

£  f*<v  -  d-o®)'1! 

k-0 
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Eq^}  -  E^Cf*)}  (9.2.9) 

and 

Eq{mJ}  -  EQ{(Y1(f*))2}  .  (9.2.10) 

For  an  element  (z,n,z* ,n')fiF,  the  value  of  the  function  g*  la  1  if  a 
passage  time  for  the  marked  job  starts  or  is  underway  when  the 
configuration  of  the  job  stack  is  z'  and  the  marked  job  is  in  position 
n* ;  the  value  of  g*  is  0  otherwise.  Formally,  let  D  be  the  state  space 
of  the  process  Z«{Z(t):t20}  appearing  in  Equation  (9.2.1).  As  in  Section  8, 
for  zeD  and  n£{l,2,. . . ,N),  we  write  h(z,n)»(i,j)  when  the  job  in  position  n 
in  the  job  stack  z  is  of  class  j  at  center  i.  Now  consider  the  embedded  jump 
chain  {V^:k20}  associated  with  the  continous  time  Markov  chain  V  of 
Equation  (9.2.2).  For  states  v'.v'^F,  the  state  spece  of  {Vfc:ki0},  we 
write  v'  r\j  v"  when  v"  is  accessible  from  v',  i.e.,  when  for  some  n£l,  the 
probability  starting  from  v'  of  entering  v"  on  the  nth  step  is  positive. 
Similarly,  for  any  subset  I  of  F  we  write  v*  <vv"  when  v"  is  accessible 
from  v'  under  the  taboo  I. 

Denoting  the  set  of  (center,  class)  pairs  in  the  network  by  C,  we  define 
a  subset  G  of  C  according  to 

G  -  {(i,j)eC:  for  some  (z,n,z' ,n')eS,  h(z* ,n')"(i,j)}  u 

{(i,j)«Cs  for  some  (z,n,z' ,n')«F-(SuT) ,  v'eS  and  v"eT  , 
v'  &  (z,n,z'  ,n’)  ,  (z,n,z'  ,n')  &v"  and  h(z'  ,n')-(i,j)}  . 


Thus,  the  set  OG^jG,,  whsre  s  (center,  class)  pair  is  in  the  set 
[respectively  G^l  if  it  is  possible  for  the  marked  job  to  be  of  class  j  at 
center  i  when  the  passage  time  specified  by  the  sets  A^,  A^,  B^,  and 
starts  [respectively  is  underway] . 

* 

Now,  for  (z,n,z* ,n')eF,  we  define  the  function  g  as 

g*(z,n,z* ,n')  -  lG(h(z',n'))  .  (9.2.11) 

(9.2.12)  PROPOSITION.  Let  g*  be  the  function  defined  by  Equation  (9.2.11), 
and  R  be  the  jump  matrix  and  £  the  vector  of  rate  parameters  for  holding 
times  in  the  continuous  time  Markov  chain  V.  Then 

E{Y1(g*)}  -  eJJ^1  g*(V(s))daj  -  (I-0R)_1(4.*oaT1)  , 

and 

E{(Y1(g*))2}  -  (1-qR)"1^* ’Sf1)  , 

where  6.  is  the  time  of  the  first  return  to  the  state  0  in  V,  and 
1 

h  -2*  »E{Y1(g  )}. 

Proposition  (9.2.12)  follows  directly  from  Theorem  (9.1.4).  We  use 
this  result  together  with  the  observation  that 
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to  obtain  the  quantities 


and 


Eo{Yi(**>} 


(9.2.14) 


EqUY^**))2}  . 


(9.2.15) 


Using  ths  ratio  formula,  Equations  (9.2.9)  and  (9.2.14)  yield  the 

2 

quantity  r.  To  obtain  the  variance  constant  a  appearing  in  the  central 
limit  theorem  (Equation  (9.2.5))  for  the  marked  job  method,  we  require  one 
more  result. 


(9.2.16)  PROPOSITION.  Let  R  be  the  jump  matrix  and  £  the  vector  of  rate 
parameters  for  holding  times  in  the  continuous  time  Harkov  chain  V.  For 
the  functions  f*  and  g*  defined  by  Equations  (9.2.7)  and  (9.2.10), 


If 


V1 

g*(V(8))d  s£  **(V 

J-0 


-  (I-nST1**  « 


*-l 


where  h'-[  (I-qR)  *£><* ")+[  (I-qR)"1^^"1)]^* -(/o^.f*) . 


Proposition  (9.2.16)  does  not  follow  directly  from  Theorem  (9.1.1)  for 
discrete  time  Markov  chains  or  from  Theorem  (9.1.4)  for  continuous  time 
Markov  chains,  but  is  established  by  similar  methods;  see  Iglehart  and 
Shadier  (1979b),  Appendix.  From  Proposition  (9.2.16)  we  obtain 
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J,(|/>i|-Eo|/061  .* 


(V(s))ds  £  f  (V  ) 
k-0 


(9.2.17) 


Then  an  expression  for  the  variance  constant  a  is 


°2  ■  pJ2l  ■  2rElli  O-J  *  ^ 


This  follows  from  Equations  (9.2.10),  (9.2.15),  and  (9.2.17). 


When  comparing  the  statistical  efflciancy  of  the  marked  job  and 
decomposition  methods,  it  is  convenient  to  have  a  central  limit  theorem 
comparable  to  Equation  (9.2.5)  but  in  terms  of  simulation  time,  t,  rather 
than  number  of  cycles,  n.  Let  m(t)  be  the  number  of  passage  times 
completed  by  time  t,  i.e.,  in  the  interval  (0,t],  If  we  denote  by  n(t)  the 
number  of  0-cycles  completed  by  time  t,  then  from  renewal  theory,  as  t-"», 

_>  —i — 
t  e^T 

with  probability  one,  where  E^Ca^}  Is  the  expected  length  of  a  0-cycle  in 
V.  This  implies  that  for  large  t,  the  number  of  0-cycles  completed  by 
time  t  is  approximately  t/Eg{ct^}.  Combining  this  result  with 
Equation  (9.2.5),  it  follows  that  as  t-*°°. 


/m(t)  \ 

■[{mOOr^E  f(Pi)j-r(f)] 
(EQ{o1})1/2a(f)/E0{M1} 


N(0,1)  . 


This  result  is  independent  of  the  initial  state  V(0).  Since  the  numerator 
in  this  central  limit  theorem  is  independent  of  the  state  0  selected  to 
form  cycles,  so  is  the  denominator.  Thus  for  the  mean  passage  time 


e  -  (E0{ot1>)1/2a/E0(M1> 

is  the  appropriate  measure  of  statistical  efficiency  for  the  marked  job 
method  and  is  Independent  of  the  state  OeS  selected  to  form  cycles.  Note 
that  we  obtain  the  quantity  Eq{ci^}  according  to 

W  •  E0{,1(1)1  > 


where  1  is  the  f  function  identically  equal  to  one  and 

Bl 

0 


E{Y 


jtt)}  -  e|J  l(V(s))ds|  -  (I-qR)"1^!*1) 


(9.2.19) 


9.3.  Variance  Constants  for  the  Decomposition  Method 

We  now  turn  to  the  decomposition  method.  As  in  Section  8.2,  we  label 
the  jobs  from  1  to  N,  and  for  i-l,2,...,N,  denote  by  Ni(t)  the  position 
of  job  1  in  the  linear  job  stack  at  time  t.  Then,  in  terms  of  the  vector 
Z(t)  corresponding  to  the  job  stack,  for  tkO  we  set 

X^t)  -  (ZCO.^Ct)) 

and 

X°(t)  -  (Z(t)  .^(t)  ,N2(t) . NN(t))  .  (9.3.1) 

Recall  that  each  of  the  processes  51«{Xi(t) :tkO)  (i-1,2, . . . ,N)  and 
X°-{x°(t):tiO}  is  an  irreducible,  positive  recurrent  continuous  time  Markov 
chain.  We  denote  the  state  space  of  the  process  £*  [respectively  X®]  by 
E1  [respectively  E0]. 
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Next  we  let  L^(t)  [respectively  L°(t)]  denote  the  last  state  visited 
by  the  Markov  chain  X  [respectively  X®]  before  jumping  to  X^(t) 

[respectively  X^(t)],  and  for  tiO  and  ial,2,...,N,  define 

Vi<t)  -  (L^O.X^t)) 

and 

V°(t)  -  (L°(t) ,X°(t) )  .  (9.3.2) 

The  process  V^“{V^(t) :t£0}  is  the  fundamental  stochastic  process  of  the 
passage  time  simulation. 

Since  each  of  the  processes  X*  and  is  an  irreducible,;  positive 
recurrent  continuous  time  Markov  chain,  so  is  each  of  the  processes 
^-{V^Ct) :t20}  and  V°.  tfe  denote  the  state  spaces  of  V1  and  V°  by  F*  and 
F^,  respectively.  The  entrances  of  V*  to  the  fixed  subset  S*  [respectively 
T*]  of  F*  (defined  by  Equation  (8.2.3))  correspond  to  the  starts 
[respectively  terminations]  of  passage  times  for  job  i.  Similarly,  the 
successive  entrances  of  to  the  fixed  subset  S®  (defined  by 
Equation  (8.2.4))  of  F^  correspond  to  the  starts  of  passage  times 
irrespective  of  job  identity,  and  the  entrances  of  V®  to  the  subset  T^ 
correspond  to  the  terminations. 

The  decomposition  method  applies  to  passage  times  for  which  the  sets  S 
and  T  (which  define  the  starts  and  terminations  of  passage  times  of  a 
particular  job)  are  disjoint.  As  in  Section  8,  {P^:n2l}  denotes  the 


128 


sequence  of  passage  times  (irrespective  of  job  Identity) ,  enumerated  in 
order  of  passage  time  start,  and  by  the  argument  in  Appendix  1,  P^—>P^. 

The  goal  of  the  simulation  is  estimation  of  the  quantity  r^(f)  of 
Equation  (8.2.5). 

Recall  that  to  estimate  r^(f),  the  decomposition  method  prescribes 

that  we  carry  out  the  simulation  of  the  process  in  random  blocks  defined 

by  the  successive  entrances  of  the  process  to  the  fixed  set  of  states  U° 

defined  by  Equation  (8.2.6).  Entrances  of  to  the  set  correspond  to 

the  terminations  of  passage  times  (irrespective  of  job  identity)  which 

occur  when  no  other  passage  times  are  underway,  and  which  leave  a  fixed 

configuration  (z^)  of  the  job  stack.  For  kkl,  denotes  the  length  In 

discrete  time  units  of  the  kth  block  (returns  to  the  set  U°)  of  the  process 

{V^:n20>;  also,  6^-0  and  6°-y!?+. .  .+y^,  mkl. 
n  U  m  i  m 

We  assume  that  V^(0)eU^  and  for  m£l  denote  by  the  number  of  passage 

m 

times  in  the  mth  block  of  the  process  V®.  Also,  we  let  Y®(f)  be  the  sum 

^  m 

of  the  quantities  f(P^)  over  the  passage  times  in  the  mth  block  of  V^, 

The  key  results  of  Section  8.2  leading  to  point  estimates  and  confidence 
Intervals  for  the  quantity  r^(f)  are  that  the  pairs  of  random  variables 

{(Y°(f),K°):m*l}  (9.3.4) 

tn  m 

are  independent  and  identically  distributed,  and  that  the  quantity 

r°(f)  -  E  {Y°(f)}/E  n{K?}  ,  (9.3.5) 

U  1  U  1 
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provided  chat  Che  quantity  E{ jf(P^) | }<«.  The  symbol  E  .{•}  is  an  abuse  of 

IT 

our  previous  notation.  It  connotes  conditional  expectation  associated 
with  starting  the  Markov  chain  in  one  of  the  states  in  the  set  U^.  The 
definition  of  the  set  implies  that  the  value  of  this  conditional 

expectation  is  independent  of  the  particular  starting  state  in  U^. 


Given  these  results,  from  a  fixed  number  of  blocks  of  V,  the 
decomposition  method  provides  the  point  estimate 

rj(f)  "  . 

n  n  n 

The  associated  confidence  interval  for  r^(f)  follows  from  the  central 
limit  theorem 


where 


n1/2{r°(f)-r°(f)> 

a°(f)/E  {K°} 

Uu  i 


->  N(0 ,1)  , 


(9.3.6) 


(a°(f))2  -  var{Y°(f)-r°(f)Kj}  . 


(9.3.7) 


Taking  f  to  be  the  identity  function,  we  restrict  attention  to  the 

quantity  r^  and  consider  computation  of  the  corresponding  variance  constant 
0  2 

(a  )  and  related  theoretical  values.  By  the  argument  which  leads  to 
Equation  (9.2.18),  the  appropriate  measure  of  the  statistical  efficiency 
of  the  simulation  is  the  quantity 

*°  “  (EuotaJ>)1/2°0'EDo^>  •  (9.3.e; 

where  a®  is  the  length  of  a  block  in  the  continuous  time  process  V°. 


130 


Th«  individual  quantities  required  to  compute  this  measure  of 
efficiency  are  defined  in  terms  of  the  successive  returns  of  the  process 
to  a  fixed  set  of  states  (U^)  rather  than  to  a  single  state.  Moreover, 
the  successive  entrances  of  to  are  not  regeneration  points  for  V^. 
Accordingly,  we  cannot  apply  the  results  of  Section  9.1  directly,  as  we 
did  for  the  marked  job  method.  Instead,  we  select  a  fixed  state 
(designated  state  u^)  from  the  set  and  compute  the  quantity  corresponding 
to  Equation  (9.3.8)  for  the  resulting  u^1- cycles .  (Note  that  the  successive 
entrances  of  the  process  to  the  fixed  state  u^  are  regeneration  points 
for  V®.)  The  expression  in  Equation  (9.3.8)  computed  for  u^-cycles  is 


e°(u°)  -  (E  Q{aJ})1/2o°/E  q{kJ}  , 
u  u 

where  the  constant  a®  (analogous  to  <3°)  is  defined  for  u°-cydes.  This 
quantity  e°(u°)  is  equal  to  e^.  To  see  this,  for  t2G  let  m°(t)  be  the 
number  of  passage  times  (irrespective  of  job  identity)  completed  in  the 
Interval  (0,t].  In  terms  of  simulation  time,  t,  we  have  the  central  limit 
theorem 


(E  {a?»1/2a°(f)/E  q{K°> 
U  u 


->  N(0,1) 


and,  when  f  is  the  identity  function,  the  variance  constant  in  the 
denominator  is  the  quantity  e*\  There  is  a  similar  central  limit  theorem 
in  terms  of  u^-cycles;  the  numerator  is  the  same  and  the  variance  constant 
in  the  denominator  is  e^(u^).  Since  the  numerators  in  these  two  central 
limit  theorems  are  the  same  as  are  the  limiting  random  variables  (N(0,1)), 


JT.  «:■!  "V  .JK* 
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tP  must  equal  e^(u^).  For  a  similar  argument,  see  Propositions  5.1  and 
5.6  of  Crane  and  Iglehart  (1975a). 

Next  we  observe  that  the  number  of  passage  times  in  a  u^-cycle  of  the 
process  V®,  as  well  as  the  sum  of  the  passage  times  in  a  u^-cycle,  does 
not  depend  on  the  identities  of  the  jobs  In  successive  configurations  of 
the  job  stack.  It  follows  that  rather  than  working  with  the  stochastic 
process  V  ,  we  can  work  with  process  W«{W(t):t20}  defined  by 

W(t)  -  (K(t).Z(t))  .  (9.3.9) 

Here  Z(t)eD  corresponds  to  the  linear  job  stack  at  time  t,  and  K(t)  is 
the  last  state  visited  by  the  Markov  chain  Z«{Z(t) :t£0)  before  jumping  to 
Z(t).  The  process  W  is  an  irreducible,  positive  recurrent  continuous  time 
Markov  chain  with  a  state  space  Chat  is  a  subset  of  D*D.  Note  that  in 
general  the  state  space  of  W  is  much  smaller  than  that  of  V®,  and  that 
working  with  the  process  W  is  computationally  advantageous. 

The  computations  rest  on  the  definition  of  two  particular  functions 

(f  and  g)  defined  on  the  state  space  of  W  and  taking  values  in  the  set 

(0,1).  We  define  the  functions  f  and  g  in  terms  of  functions  f®  and  g® 

defined  on  F^,  the  state  space  of  the  process  V^.  We  take  the  function 

f®  to  be  the  indicator  function,  1  n,  of  the  set  which  defines  the 

SU 

starts  of  passage  times  irrespective  of  job  identity,  i.e.,  for 
(z,nl . nN,z',ni’ . nN)  €F°  » 

f  (z »o^»  •  •  •  •  '  «ni»  •  •  •  *°ij)  "  1  q(*  *n2_,  •  •  •  »°jj)  •  (9.3.10) 
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Thus  if  a  passaga  time  for  some  job  starts  when  hits 

(z.n^ . r^.z' ,n|, . . . ,n^) ,  then  f°«l.  Note  that  for  each  (z,z')  in  the 

state  space  of  W,  there  exist  n^,...,n^,n^ .  and  n^  such  that 

(z.n^, . . . »nN,z* ,n| . n^)€F°.  For  a  state  (z,z')  of  W,  we  define 

f (z » z 1 )  -  f°(z,n1,...,nN,z' ,n^,...,n^)  .  (9.3.11) 

The  function  f  is  well  defined  since,  for  fixed  z  and  z',  the  function  f 
is  independent  of  its  other  arguments. 


For  an  element  (z.r^ . n^.z* ,n^, . . . ,n^)eF°,  the  value  of  the  function 

is  the  number  of  passage  times  that  start  or  are  underway  when  the 
configuration  of  the  job  stack  is  z'.  Formally,  for 
(z,n1,...,nN,z,,n;[,...,n^)€F0,  we  define 

N 


g  (z,nx . njj.z  ,^,...,11^)  m  2J  lg(h(z,k))  . 

Then,  for  (z,z')  in  the  state  space  of  W,  we  define 


(9.3.12) 


g(z,z') 


0/ 

g  (*.«-, 


’V 


•«i) 


(9.3.13) 


The  justification  for  using  the  process  W  is  that  the  number  of  passage 
times  (which  start  and  terminate)  in  the  first  u^-cycle  of  is 


V1 


£  f(w  )  , 


3-0 


(9.3.14) 


and  the  sum  of  the  passage  times  in  the  first  u^-cycle  of  is 


g(W(s))ds  . 


(9.3.15) 
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Here  (respectively  K^)  Is  the  time  of  the  first  return  of  the  process  U 
(respectively  the  jump  chain  {W^:k£0})  to  the  fixed  state  w®.  The  return 
state  w®  corresponds  to  the  fixed  state  u^  selected  from  the  set  U®,  i.e., 
if 


"0  *  ‘‘•“I . °N’*0,ni . "i> 


then  w^-(z,z®). 


Proposition  (9.3.16)  follows  directly  from  Theorem  (9.1.1). 

(9.3.16)  PROPOSITION.  Let  f  be  the  function  defined  by  Equation  (9.3.11) , 
and  R  be  the  transition  matrix  of  the  discrete  time  Markov  chain  {W.  :k^0). 
Then 

1 

and 

E{(Yx(f))2}  -  (I-qR)'1^  , 

where  is  the  time  of  the  first  return  to  the  state  w°  in  {Wfc:k£0}  and 
h-2f oE{YL(f ) }-f 

From  Proposition  (9.3.16)  we  obtain  the  quantities  E  q{K^}  and  E  ^{(K^)2} 

u  u 

according  to 

E  0{K1}  “  E  0{Yl(f>}  (9.3.17) 

u  w  A 


E(Y1(f)}  -  Ei 
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and 

E  Q{ (K?)2}  -  E  Q{(Y  (f))2}  .  (9.3.18) 

u  wu 

(9.3.19)  PROPOSITION.  Let  g  be  Che  £uncClon  defined  by  Equation  (9.3.12), 
and  R  be  the  jump  matrix  and  £,  the  vector  of  rate  parameters  for  holding 
times  in  the  continuous  time  Markov  chain  W.  Then 

E(Y1(g)}  -  e|  J^1  g(W(s))ds|  -  (I-qR)"1^’1)  . 

and 

E((Y1(g))2}  -  (I-qR)-1^'1)  , 

where  ^  is  the  time  of  the  first  return  to  the  state  w°  in  W,  and 
h-2f°E{Y1(g)}. 


Proposition  (9.3.19)  follows  directly  from  Theorem  (9.1.4).  Ue  use 
this  result  to  obtain  the  quantities 


E 

u 


E  p 


°(j=i 


E  0CY1Cg)> 
w 


(9.3.20) 


and 


E  Q{(Y1(g))2}  . 


(9.3.21) 


Using 

Analogous 


the  ratio  formula.  Equations  (9.3.17)  and  (9.3.20)  yield  r®. 
to  Proposition  (9.2.16)  we  have 


_  .  rinoiii — -tfnrT"TTrnmrgiriiMm 
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(9.3.22)  PROPOSITION.  For  EtY^f)}  and  Ett^g)}  given  by 
Propositions  (9.3.16)  and  (9.3.19), 


jj^1  g(W(s))ds  22  f(Wk)|  -  (I-qR)"1!!  , 


where  h-(A«£_1)  ®E{Y1(f )  }+f  »E{Y1(g)  of . 


We  use  this  result  to  obtain 


BoKe'pJWI-eJ/51  g(W(s))ds  2  3 
UU(\j-0  V  WU(J0  k-0 


f(wk)  . 


(9.3.23) 


Then  to  compute  the  variance  constant  (a  )  ,  we  use  the  results  o£ 
Equations  (9.3.18),  (9.3.21),  and  (9.3.23). 


9.4.  Numerical  Results 

We  once  again  consider  the  closed  network  of  queues  of  Section  5.1, 
and  the  limiting  passage  times  P  and  R  therein.  Recall  that  the  limiting 
passage  time  P  starts  when  a  job  joins  the  center  1  queue  upon  completion 
of  a  center  2  service  and  terminates  when  the  job  next  joins  the  center  2 
queue.  Similarly,  the  response  time  R  is  associated  with  the  time  between 
successive  entrances  of  a  job  into  the  center  1  queue  upon  completion  of 
a  center  2  service. 


For  the  passage  time  P,  the  sets  A1  and  A 2  defining  the  starts  of 
passage  times  are 
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^  -  { (i,N) :0S1<N}  , 

and 

A2  -  { (1,1) :0<iSN}  . 

Similarly,  the  sets  and  B2  defining  the  terminations  of  the  passage 
time  P  are 

B1  -  { (1,1) :0<l£N} 

and 

B2  -  { (1-1,1) :0<i£N}  . 

For  the  response  time  R,  the  sets  A^  and  A2  are  the  same  as  for  the  passage 
time  P,  but  B^-A^  and  B2"A2‘ 

In  connection  with  the  marked  job  method,  the  process 
V"{ (L(t) ,X(t)) :t20),  where  L(t)  is  the  last  state  visited  by  the  Markov 
chain  X  before  jumping  to  X(t),  has  state  space 

F  -  { (i,j ,1+1, j+1) :0Sj<N,  lsi<N)  u  { (i,N, 1+1,1) :0Si<N>  u 

{(i,j  ,1-1, J)  :0<i£N,  lSjSN}  u  {(1,1, 1,1)  :K1SN}  . 

The  subsets  of  F  defining  the  starts  and  terminations  of  passage  times 
for  the  marked  job  are 


S  -  {(i,N, 1+1,1) :0S1<N> 
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and 


T  -  {(1,1, 1-1,1) :0<iSN}  . 

Tables  9.1  and  9.2  give  theoretical  values  for  simulation  of  the  closed 

network  of  queues  by  the  marked  job  method.  Numerical  results  are 

displayed  for  the  mean  of  the  response  time  R  and  corresponding  results 

for  the  passage  time  F  are  in  parentheses.  For  the  case  of  N-2  jobs 

(Table  9.1),  the  set  S-{0,2,1,1),  (1,2, 1,1)}.  With  Aj-1,  X2«-0.5,  and  p-0.75, 

the  numerical  results  show  that  on  the  average  O-cycles  defined  by  returns 

to  the  state  (0,2, 1,1)  are  twice  as  long  as  those  defined  by  the  returns 

2 

to  the  state  (1,2, 1,1).  Note  that  as  expected,  the  quantities  a  /Eq{M^} 

1/2 

(as  well  as  e«(Eg{a^})  o/Eq{M^})  are  the  same  for  the  two  return  states. 

Table  9.2  gives  results  for  N»4  jobs.  Here  there  are  four  possible  return 
states,  and  for  the  parameter  values  selected,  returns  to  the  state 
(3, 4, 4,1)  occur  most  frequently,  and  on  the  average  eight  times  more  often 
than  returns  to  the  state  (0,4, 1,1). 

We  now  turn  to  the  decomposition  method.  As  we  saw  in  Section  8.2, 
the  process 

X°  -  {(Z(t),N1(t),...,NN(t)):t£0} 

has  state  space  E®,  where 

E°  -  {(i,^, . . .  ,njj)  :0SiSN;  ls^, . . .  .n^N;  n^j  for  ii*j}  . 

The  underlying  continuous  time  process  V®  defined  by 
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V°(t)  -  (L°(t),X°(t))  , 

where  L^(t)  is  the  last  state  visited  by  the  Markov  chain  J®  before  jumping 
to  X^(t) ,  has  state  space  F®.  The  subsets  of  F®  defining  starts  and 
terminations  of  passage  times  are 

S®  ■  { (i,n^, . . . .n^.i+l.n^, . . . .n^) :0£i<N;  for  exactly  one  j*, 

n  #-N  and  n'*»l;  lstaj<N  end  nj-n^+1,  j^j*;  n^J^  for  Wj} 

and 

T  ■  { (i ,n^  , • • • >n^, 1~1 ,n^ j « ■ i b^) i 0^l£N  j 

lSn^-nj^N,  for  1<: j 5N ; 

The  process  Z*{Z(t) :tiO}  has  state  space  D"{0,1, . . . ,N} ,  the  set  D°«{0}, 
and  the  set  17®  defining  blocks  of  the  process  V®  is 

U  *  {  (ltn^i  •  •  •  tn^iOiD^!  •  >  ■  |Bjj)  !l£n^«  •  >  •  >n^^Ki  n^^n^  for  ij^j}  > 

The  state  space  of  the  stochastic  process  W  is 

{(i,i+l)  :0S1SH-1}  u  { (i,i-l)  :lSis;N}  , 

and  the  state  w®»(l,0). 

Table  9.3  gives  theoretical  values  for  simulation  of  the  closed  network 
of  queues  by  the  decomposition  method  for  the  mean  of  the  passage  time  P. 
The  table  gives  results  for  N*1  to  N»4  jobs,  and  the  parameter  values  are 
the  same  as  in  Tables  9.1  and  9.2.  For  N-2  jobs,  the  value  of  the 
quantity  e®»e®(u®)  of  Equation  (9.3.8)  which  measures  the  statistical 
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efficiency  of  the  decomposition  method  is  16.546.  The  corresponding  value 
from  Table  9.1  for  the  marked  job  method  is  20.890.  Thus,  for  these 
parameter  values  the  decomposition  method  is  approximately  21  percent  more 
efficient  than  the  marked  job  method.  For  N“4  jobs,  the  decomposition 
method  is  41  percent  more  efficient. 


Numerical  results  bearing  on  the  statistical  efficiency  of  the 
decomposition  method  for  simulation  of  the  closed  network  of  queues  appear 
in  Table  9.4.  For  N>1  to  Na6  jobs,  the  table  gives  theoretical  values  of 
the  quantities  and  e^  for  three  sets  of  parameters  values.  We  hold  the 
value  of  ^2*1  and  p*0.75  fixed,  but  vary  X^.  Table  9.5  gives  a  comparison 
of  the  relative  efficiency  (e/e®)  of  the  marked  job  and  decomposition 
methods  for  the  same  sets  of  parameter  values. 
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TABLE  9.2 

Theoretical  Values  for  the  Marked  Job  Method. 
Passage  Time  R  (P)  In  Closed  Network  of  Queues. 
N«4,  X^»1.0,  p«0.75. 


Parameter 

Return  State  of 

< 

ft 

w 

ft 

IV 

o 

(0,4, 1,1) 

(1,4, 2,1) 

(2,4, 3,1) 

(3, 4, 4,1) 

W 

216.0 

108.0 

54.0 

27.0 

(M1  ) 

E°l?l  PJ) 

248.0 

(196.0) 

124.0 

(98.0) 

62.0 

(49.0) 

61.0 

(24.5) 

W 

15.0 

(15.0) 

7.5 

(7.5) 

3.75 

(3.75) 

1.875 

(1.875) 

ft  ) 

eo{2  pjj/Eoft> 

16.533 

(13.067) 

16.533 

(13.067) 

16.533 

(13.067) 

16.533 

(13.067) 

o2 

2111.343 

(2139.600) 

1055.672 

(1069.800) 

527.836 

(534.900) 

263.918 

(267.450) 

a2/EQ{M1} 

140.756 

(142.640) 

140.756 

(142.640) 

140.756 

(142.640) 

140.756 

(142.640) 

('oft1) 1/2  o/E0<Ml> 

48.241 

(48.562) 

48.241 

(48.562) 

48.241 

(48.562) 

48.241 

(48.562) 
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TABLE  9.3 

Theoretical  Values  for  the  Decomposition  Method. 
Passage  Time  P  in  Closed  Network  of  Queues. 
Xj*1.0,  Xj^O.S,  p«0.75. 


Parameter 


E  01“°1> 
U 


«o  E  p? 

u  (j"i  J 


196.0 


Eo<*?> 

U 


Eu°(j5i  PJ)/Eu°{Kl 


6.667 


9.714  13.067 


,-0.2 


176.0  1023.673  4317.227 


CaJ)2/E  0{kJ} 

u 


58.667  146.249  287.815 


(E  olai})1/2  °o/E  olEi> 

\  U  /  U 


9.798  16.546  25.035 


34.491 


10.0.  NETWORKS  WITH  MULTIPLE  JOB  TYPES 

We  have  considered  in  previous  sections  the  problem  of  simulating 
closed  and  finite  capacity  open  networks  of  queues,  respectively,  for 
general  characteristics  of  passage  times.  Under  consideration  here  are 
networks  with  multiple  job  types  and  the  estimation  of  individual  and 
joint  characteristics  of  passage  times  over  the  several  job  types.  The 
type  of  a  job  may  influence  its  routing  through  the  network  as  well  as 
its  service  requirements  at  each  center.  For  expository  convenience,  we 
assume  that  there  are  only  two  job  types  in  the  network  and  we  mark  one 
job  of  each  type.  By  tracking  these  two  jobs,  we  are  able  to  produce  from 
a  single  replication  confidence  intervals  for  a  variety  of  passage  time 
characteristics.  The  estimation  method  of  this  section  can  also  be  applied 
to  networks  with  only  a  single  job  type;  the  result  is  an  alternative 
scheme  to  that  proposed  in  Section  7. 

10.1.  Preliminaries 

We  consider  closed  networks  of  queues  with  a  finite  number  of  jobs, 

N,  of  two  types,  and  assume  that  there  are  [respectively  N^]  jobs  of 
type  1  [respectively  type  2]  with  N^+N^"**.  In  each  network  there  are  a 
finite  number  of  service  centers,  s,  and  a  finite  number  of  job  classes, 
c.  All  jobs  retain  their  job  type,  but  may  change  class  as  they  traverse 
the  network.  (Think  of  type  1  jobs  as  cubes  and  type  2  jobs  as  spheres, 
and  let  job  classes  correspond  to  different  colors.  Then  we  permit  jobs 
to  change  color,  but  not  shape.)  Upon  completion  of  service  at  center  i, 
a  type  v  job  of  class  j  goes  to  center  k  and  changes  to  class  i,  with 
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probability  We  assume  that  for  v»l,2,  £^V^*{p^k^:lsi,kis,  lsj,Jl5c} 

is  a  given  irreducible  Markov  matrix. 

The  service  times  and  service  discipline  at  each  service  center  are 
as  in  Section  7  with  the  exception  that  they  may  also  depend  on  job  type. 

We  briefly  review  the  situation.  At  each  service  center  jobs  queue  and 
receive  service  according  to  a  fixed  priority  scheme  among  classes  and 
types,  which  scheme  can  vary  from  center  to  center.  Each  center  operates 
as  a  single  server,  processing  jobs  of  a  fixed  type  and  class  according 
to  a  fixed  service  discipline.  All  service  times  in  the  network  are 
mutually  independent,  and  at  each  center  have  a  distribution  with  a 

Cox-phase  representation  with  parameters  which  may  depend  on  the 
service  center,  type  and  class  of  job  being  serviced,  and  the  "state"  of 
the  entire  system.  (As  usual,  we  exclude  zero  service  times  occurring 
with  positive  probability.)  A  job  in  service  may  or  may  not  be  preempted 
(according  to  a  fixed  procedure  for  each  center)  if  another  job  of  higher 
priority  joins  the  queue  at  the  center. 

We  restrict  the  present  discussion  to  networks  in  which  all  service 
times  are  exponentially  distributed,  and  deal  with  distributions  having 
a  Cox-phase  representation  in  the  usual  way  by  the  method  of 
stages.  To  characterize  the  state  of  the  system  at  time  t,  we  let  S^(t) 
denote  the  (type,  class)  pair  of  the  job  receiving  service  at  center  i  at 
time  t,  i*l,2,...,s.  If  there  are  no  jobs  at  center  i  at  time  t,  we  set 
Si(t)-(0,0).  We  denote  by  j^i), . . .  .j^^  (i)  the  (type,  class)  pairs 
served  at  center  i  ordered  by  decreasing  priority,  and  let 


147 


cJ4)Ct) . cj1>  (t)  denote  the  number  of  jobs  in  queue  at  time  t  of  the 

J1  Jk(i) 

various  (type,  class)  pairs  served  at  center  i.  We  mark  one  job  of  each 
of  the  two  types  in  order  to  measure  their  passage  or  response  times.  As 
in  previous  sections,  we  view  the  N  jobs  as  being  completely  ordered  in  a 
linear  stack,  and  let  the  vector  Z (t)  be  given  by: 


Z(t)  -  (C.(1)  (t),...,C,(1)(t),S.(t);...; 

Jk(l)  h  L 


C.Cs)  (t),...,C,(s)(t),S  (t))  . 
jk(s)  jl  8 


(10.1.1) 


The  linear  job  stack  again  corresponds  to  the  order  of  components  in  the 
vector  Z(t)  after  ignoring  any -zero  components.  Within  a  (type,  class) 
pair  at  a  center,  jobs  waiting  appear  in  the  job  stack  in  the  order  of  their 
arrival  in  the  center,  the  latest  to  arrive  being  closest  to  the  top  of 
the  stack.  Let  N^(t)  (v-1,2)  denote  the  position  from  the  top  in  this 
job  stack  of  the  type  v  marked  job.  Then  for  tkO,  the  state  vector  of 
the  network  is 


X (t)  -  (Z(t),N1(t),N2(t)) 


(10.1.2) 


Under  the  exponential  service  time  and  Markovian  routing  assumptions,  the 
process  X*{X(t) :t<:0}  is  an  irreducible  continuous  time  Markov  chain  with 
finite  state  space  E. 


10.2.  Simulation  for  Passage  Times 


We  specify  the  passage  (or  response)  times  for  the  two  types  of  jobs 
by  eight  subsets  of  E:  for  v-1,2.  The  sets 

and  [respectively  B^ ,  B^v*]  determine  when  to  start  [respectively 


fc 


HUMS.  -  « 


stop]  the  clock  measuring  a  particular  passage  time  for  the  type  v  marked 
job.  Denoting  the  jump  times  of  the  process  X  by  {TQ:nkO},  for  k,n2l,  we 
require  that  the  sets  A^ ,  A^ »  and  satisfy 

if  X(Tn-1)eAjv),  X(Ta)£A^V),  X(TQ_1+k)eA^V)  and  X(TQ+k)€A<V)  , 

then  X(T  . .  )€B.[V^  and  X(t _ )gBjV^  for  some  0<m£k  ; 

o—iTtn  i.  nna  i 


if  X(Vl).B<V>,  X(Tn).B‘V>,  X<V1+k).B'V>  *nd  X(W).B<V)  , 

then  X(T  , .  )«A^V^  and  X(T  .  )«a5V^  for  some  0Sm<k  . 
n-l-Hn  l  tiTtn  z 

Also,  in  terms  of  the  jump  times  of  X,  we  define  four  sequences  of  random 
times:  {sjV^:jS0}  and  {TjV^:j£l},  for  v«l,2.  The  start  [respectively 
termination]  time  of  the  jth  passage  time  for  the  type  v  marked  job  is 
denoted  by  [respectively  TjV^].  Formally,  we  have  for  v-1,2, 

sj'0  -  inf (VTj(V)  ;X(Tn)«A^V)  ,  XtT^JeA^},  j*0 

T<V)  -  inf{Tn>Sj(^)1:X(Tn)eB^V),  XCT^eB*'0 } ,  j*l  . 

(v)  (v)  (v) 

The  jth  passage  time  for  the  type  v  marked  job  is  -T^  '-Sj_£,  jkl. 
Note  that  the  definition  of  these  times  is  as  in  Section  7.  For  response 
times  of  type  v  jobs,  A^*B^V\  A^^B^,  and  sjV^*TjV^  for  all  J£l. 


Let  L(t)  denote  the  last  state  visited  by  the  Harkov  chain  £  before 
jumping  to  X(t) ,  and  for  t*0  set 


V(t)  -  (L(t) ,X(t))  . 


(10.2.1) 
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The  process  V*(V(t):t£0}  has  a  state  space  F  consisting  of  all  pairs  of 
states  (i,j),  i,j€E,  for  which  a  transition  in  X  from  state  i  to  state  j 
can  occur  with  positive  probability.  In  general,  of  course,  the  size  of 
the  state  space  F  is  larger  than  that  of  E.  The  "Q-matrix"  used  in 
generating  the  Markov  chain  V  can  be  obtained  easily  from  that  for  £. 
Since  J  is  an  irreducible,  positive  recurrent  Markov  chain,  so  is  V. 
Clearly,  the  entrance  times  of  V  to  a  state  (i,j)«F  correspond  to  the 
times  of  transition  in  X  from  state  1  to  state  j.  For  a  type  v  job,  we 
define  two  subsets  of  F  according  to: 

S(V)  -  {(i,j)eF:ieA<V) ,  jeA<v)} 

T(V)  -  {(i,j)£F:i£BjV>,  j£B^v)}  . 

Thus  the  entrances  of  V  to  S  ^  [respectively  T^J  correspond  to  the 
start,  [respectively  termination]  times  of  passage  times  for  the  type  v 
marked  job.  Of  course  for  response  times  of  a  type  v  job,  S^-T^. 

The  argument  employed  in  Appendix  1  shows  that  for  v»l,2,  the 

(v)  (v) 

sequence  converges  in  distribution  to  a  random  variable  P 

Moreover,  the  sequence  of  passage  times  of  type  v  jobs  (irrespective  of 

job  Identity)  in  the  order  of  start  (or  termination)  also  converges  in 

distribution  to  P^v\  Our  concern  is  with  the  estimation  of 

characteristics  associated  with  these  limiting  passage  times. 
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Estimation  of  E{R^}  and  P{R^£x} 

Using  the  process  V  defined  by  Equations  (10.1.1),  (10.1.2)  and 
(10.2.1),  we  consider  first  the  estimation  of  characteristics  of  the 
limiting  response  time,  R^,  of  a  type  1  job.  For  this  estimation 
problem,  of  course,  it  is  not  necessary  to  mark  a  type  2  job.  Since  R^ 
is  a  response  time,  We  select  a  fixed  state  of  S^,  which 

for  convenience  we  designate  state  0,  and  assume  that  V(0)«0. 


Suppose  that  we  wish  to  estimate  E{R^  '}.  The  successive  entrances 

of  V  to  S^  constitute  the  starts  and  terminations  of  response  times  of 

the  type  1  marked  job.  Let  R^  (n20)  denote  the  time  between  the  nth 

and  (n+l)st  entrances  to  S^\  with  the  0th  entrance  to  S^  occurring  at 

t-0.  Also,  let  {VQ:n20}  denote  the  embedded  jump  chain  associated  with 

V.  The  random  times  {<XQ:n£l}  and  {y^nkl}  denote  the  lengths  of  the 

successive  0-cydes  (successive  returns  to  the  fixed  state  0)  for  V  and 

{V  :n20} ,  respectively.  Then  the  number  of  response  times  for  the  type  1 
n 

marked  job  in  the  first  0-cycle  of  V  is 

"i”  '  i §> 

where  6«“0  and  6  *y,+...+Y  >  mil.  The  sum  of  the  response  times  in  that 
u  ml  m 

cycle  is  simply 


a  rn  T  rW 

1  n-1  n  ' 

We  denote  the  analogous  quantities  in  the  kth  0-cycle  by  and  a^.  The 
fact  that  Visa  regenerative  process ,  together  with  a  renewal  argument 
(cf.  Appendix  2)  establishes 
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(10.2.2)  PROPOSITION.  The  pairs  of  random  variablas  {(ot^N^)  :lul)  ara 
independent  and  Identically  distributed.  Provided  that  E{R^}<“, 

E{R(1)}  -  E{a1>/E{Nj1)}  . 


At  this  point,  the  arguments  of  the  standard  regenerative  method 

—  — (1) 

hold  and,  based  on  n  cycles,  we  can  construct  the  point  estimate  a^/lT  ' 

2 

and  (provided  that  an  estimate  is  available  for  a  ,  the  variance  of 
aj-E{R^}N^)  an  associated  confidence  interval  for  E{R^}.  The 
confidence  Interval  is  obtained  from  the  central  limit  theorem 


n1/2ta  /N(1)-E(R(1>>] 
n  n 

a/EtN^h 


N(0,1) 


Here  a  -(a.,+.  ..+a  )/n  and  N^1^-(N.^1V. .  .+N^)/n. 
n  l  n  n  i  n 


If  we  are  interested  in  the  distribution  function,  P{R^£x}  of  R^ , 
we  proceed  as  above,  but  define  in  addition  the  1.1. d.  sequence  of  random 
variables  {Y^:k2l},  where,  for  example, 

N<X> 

Yl  '  ti  ' 

Then  the  point  estimate  of  P{R^£x)  is  just  Yq/N^  ,  and  we  obtain 
confidence  intervals  in  the  usual  way. 


Estimation  of  E{R^}  and  E{R^} 

Now  suppose  that  we  wish  to  estimate  the  expected  passage  time  for 
type  2  Jobs,  E{R^},  as  well  as  E{R^}.  Response  times  for  the  type 
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marked  job  start  and  terminate  at  the  entrance  times  of  V  to  the  set 
S^-T^.  Let  N^2^  denote  the  number  of  entrances  to  S^  of  V  in  the 
kth  0-cycle.  For  example,  in  the  first  0-cyde 

«rl 

<2)  ■  E  1{v  eS«>>  • 

n-0  n 

Although  we  are  able  to  begin  the  simulation  at  the  start  of  a  response 

time  for  the  type  1  marked  job,  in  general  a  response  time  for  the  type  2 

marked  job  is  underway  at  time  t«0.  Similarly,  at  the  end  of  a  0-cyde, 

a  response  time  for  the  type  1  marked  job  terminates,  but  a  response  time 

for  the  type  2  marked  job  is  still  underway.  After  n  0-cycles, 

(2)  ( 2 ) 

Nj  response  times  for  the  type  2  marked  job  have  started  and 

the  sum  of  these  response  times  is  approximately  0^+. . •+ctn*  The  error  in 

this  approximation  is  due  to  the  partial  response  time  at  t»0  which  is 
(2)  (2) 

not  counted  in  N,  . .+NV  and  the  last  response  time  which  is  counted, 
x  n 

but  does  not  terminate  before  the  end  of  the  nth  0-cycle.  Since  the  point 
estimates  and  confidence  intervals  here  are  based  on  large  sample  theory 
(strong  laws  and  central  limit  theorems) ,  these  errors  are  negligible  for 
n  large.  In  fact,  the  errors  due  to  the  two  response  times  at  t-0  and  at 
the  end  of  the  simulation  run  compensate  for  each  other.  Consequently, 
we  have 


(10.2.3)  PROPOSITION.  The  pairs  of  random  variables  {(ci^,N^):k2:l}  are 

(2) 

independent  and  identically  distributed.  Provided  that  E{RV  }<°°, 

E{R<2)}  -  E{a1}/E{N^2)}  . 
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In  the  presence  of  Proposition  (10.2.3),  the  point  estimate  of  E{R^}  is 
—  —(2) 

a  /N  ,  end  we  can  use  the  standard  regenerative  method  to  obtain  a 
n  xi 

confidence  interval. 


Estimation  of  E{R(1)}-e{R(2) } 

Suppose  now  that  we  wish  to  estimate  r^-r^2\  where  r^«E{R^} 

and  r'  '-E(RV  We  can  take  as  a  point  estimate  the  quantity 

(a  /N^)-(a  /N^),  but  need  a  bivariate  central  limit  theorem  in  order  to 
n  n  n  n 

produce  a  confidence  interval.  To  this  end,  we  let 


and 


for  kkl.  (We  take  all  our  vectors  to  be  column  vectors.)  The  random 
vectors  {Z^:k2l}  are  l.i.d.  since  each  Z^  is  only  a  function  of  the  kth 
0-cycle.  Furthermore,  Equations  (10.1.4)  and  (10.1.5)  imply  that  E{Z^}-£. 
Denoting  the  transpose  of  by  Z£,  let  Z"E{Z^Z^}-{o^^ }  be  the  covariance 
matrix  of  the  Z^' s .  Assuming  that  the  elements  of  £  are  finite,  we  have 
the  central  limit  theorem 


n-1/2  E  Z.^  N(0,Z)  ,  (10.2.4) 

k-1 

where  N(£,Z)  is  a  multivariate  normal  random  variable  with  zero  mean  vector 
and  covariance  matrix  Z.  We  can  rewrite  Equation  (10.2.4)  in  the  form 


(N<2) /E{n{1) }) E{N(1) } [ { (T  /W*15 )-r (1) } ] 

Q  1  1  Q  u 

(N<2)/E(N^2)})E{nP)>[{(o  /i?J2))-r(2)}} 
u  i  ±  an 


->  N(0,Z)  . 


(10.2.5) 
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Since  E{Nf^}/N^->l,  ve  can  use  Lemma  (2.1.8)  Co  conclude  that  Che 
x  n 

factors  (N<v)/E{N^})  in  Equation  (10.2.5)  can  be  dropped.  With  these 
factors  removed,  again  apply  Lemma  (2.1.8)  with  the  mapping  h  given  by 


h(xrx2)  -  (x1/E{N^1)>,  x2/E{n[2)}) 


to  obtain 


n1/2 


)-r 

)-r 


(1) 

(2) 


— > 


N(O.BZB')  , 


(10.2.6) 


where 

l/EtN^15}  0 

I  -  1  fn)  • 

o  i/e{nJ2)} 

Note  that  from  Equation  (10.2.6)  we  could  construct  a  simultaneous 
confidence  interval  for  (r^  ,r^) .  Finally,  a  third  application  of 
Lemma  (2.1.8),  this  time  using  h(x^,x2)*x^-x2,  yields 


(10.2.7)  PROPOSITION.  Provided  that  a22<00, 

(n1/2/a)[{(a_ /N^1))-(a V2))}-(r(1)-r(2))]  ->  N(0,l)  (10.2.8) 

u  n  n  n 

where 

°2  “  o11/*2H'i1>>  +  o22/e2{n£2)>  -  2o12/(e{nJ1)}e(n^2)})  . 

We  can  use  the  central  limit  theorem  of  Equation  (10.2.8)  to  construct  a 
confidence  interval  for  r^-r^2\  provided  that  an  estimate  for  the 
constant  a  is  available.  Using  the  classical  method,  we  can  estimate  o 
from  the  sequence  of  observations  taken  in  the  n  0-cycles  of  the  process  V. 
This  estimate  for  a  appears  in  Appendix  3. 
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ff 

Since  6^(1, 0)  is  one  possible  value  of  j5,  using  £  is  guaranteed  to  yield 

a  variance  reduction  over  that  obtained  by  marking  just  one  job.  Again, 

2  * 

of  course,  we  must  estimate  the  variance  cr(£  )  given  in  Equation  (10.2.9) 
from  the  observations  recorded. 


Estimation  of  P{R(1)Sx}-P{R(2)sx} 

Finally,  we  consider  the  estimation  of  P{R^£x}-P{R^£x)  for  a  given 
value  of  x.  This  is  the  most  difficult  of  the  problems  for  networks  with 
multiple  job  types  that  we  treat.  Since  the  value  of  x  is  fixed  throughout 
the  discussion,  in  general  we  suppress  in  our  notation  the  dependence  of 
x.  Again  we  form  0-cycles  based  on  the  response  times  for  the  type  1 
marked  job.  Here,  however,  when  a  0-cycle  ends,  we  do  not  know  whether 
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Che  response  time  for  Che  type  2  marked  job  in  progress  will  be  less  Chan 
or  equal  Co  x.  Thus,  with  respect  Co  the  response  times  for  the  type  2 
marked  job,  the  0-cycles  used  previously  do  not  create  the  i.i.d.  cycles 
needed  Co  establish  a  central  limit  theorem.  Instead,  we  form  new  cycles 
by  grouping  together  a  random  number  of  consecutive  0-cycles.  Let 
tj*cij+.  .  .a^,  ial.  Then  let  Sj^  be  the  start  time  of  the  response  time  for 
the  type  2  marked  job  underway  at  the  conclusion  of  the  ith  0-cycle.  We 
assume  that  V(0)*0  and  regard  the  value  of  the  response  time  for  the  type  2 
marked  job  underway  at  the  start  of  the  simulation  to  be  greater  than  (the 
fixed)  x.  We  do  this  so  that  the  start  of  the  simulation  corresponds  to 
the  beginning  of  one  of  the  new  "super-cycles"  we  are  constructing. 

Defining  a  random  variable  y  according  to 


Y  “  inf {i2l:ti~si>x}  , 

the  length  of  the  first  super-cycle  is  simply  a^+a^. .  ,+a^,  and  the  number 
of  response  times  for  the  type  v  marked  job  started  in  this  super-cycle 
is  n^V^»N^V^+. . .+N^V\  Successive  super-cycles  are  defined  in  an  analogous 
fashion.  For  k2l,  we  define  to  be  the  number  of  response  times 

terminating  in  the  kth  super-cycle  which  are  less  than  or  equal  to  x; 


°1  1 


Observe  that  by  the  definition  of  a  super-cycle,  the  first  response  time 
of  the  type  2  marked  job  terminating  within  a  super-cycle  must  be  greater 


than  x.  Thus  we  have 
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(10.2.10)  PROPOSITION.  The  random  variables  (Y^2^  :ksl)  are  independent 
and  identically  distributed. 


Of  course,  the  Y^^'s  are  i.i.d.  also,  as  are  the  n^^'s  and  n£2^'s.  We 
can  now  form  the  bivariate  central  limit  theorem  analogous  to  Equation 
(10.2.6),  namely 


1/2 


(Y^15  /n^  )-P{R(1)  Sx} 

(7(2)/^2)).p{R(2)Sx} 


->  N(0,B(x)£(x)B' (x))  , 


where 


B(x)  - 


l/Etn^} 


1/E(n^2>> 


and 


with 


Z(x)  -  {a^/x)} 


(x)  -  E{[Yji)-n^i)P{R(i)Sx}I[Y^)-n^)P{R(J^x}]}  . 


Finally,  by  the  same  argument  used  in  Proposition  (10.2.7),  we  obtain 


(10.2.11)  PROPOSITION.  Provided  that  au(x),  o12(x),  a22(x)<«, 

(n1/2/o(x))[(Y(1)/n(1)-Y(2)/n(2))-(P{R(1)sx}-P{R(2)Sx})]  ->  N(0,1)  ,  (10.2.12) 
n  n  n 

where 

o2(x)  -  0n(x)/E2  n£1J}  +  o22(x)/E2{n^2)}  -  2o12(x)/(E{n<1) }E{n<2)})  . 


We  can  estimate  the  quantity  a(x)  from  the  observations  in  the 
n  super-cycles  using  the  classical  method;  see  Appendix  3.  Then  we  construct 
confidence  intervals  for  P{R^£x}-P{R^Sx}  from  Equation  (10.2.12)  in 
the  usual  way. 


The  discussion  in  this  section  has  concentratad  on  problems  associated 
with  the  estimation  of  characteristics  of  response  times  for  the  two  types 
of  jobs.  The  estimation  of  characteristics  of  two  passage  times,  or  one 
response  time  and  one  passage  time,  is  in  general  easier.  This  is  because 
there  is  the  possibility  of  forming,  from  0-cycles  based  on  one  type  of 
job,  super-cycles  which  terminate  when  no  passage  time  of  the  other  type 
of  job  is  underway. 


We  have  considered  explicitly  only  the  case  of  two  job  types.  The 
estimation  methods  of  this  section  apply  equally  well  to  networks  having 
more  than  two  job  types.  The  state  space  which  results  from  the 
augmentation  of  the  vector  X(t)  (by  components  to  track  a  marked  job  of 
each  of  the  job  types)  is  of  course  larger. 


10.3.  Example  and  Numerical  Results 

To  illustrate  the  technique  of  the  previous  section  for  estimation  of 
response  times,  we  consider  a  simple  closed  network  of  queues  having  two 
types  of  jobs  and  two  service  centers;  see  Figure  10.1.  There  are  N  jobs 
in  the  network,  N^  jobs  of  type  1  and  N^  jobs  of  type  2.  After  completion 
of  service  in  center  1,  a  type  v  job  joins  the  queue  at  center  1  (with 
probability  p^)  or  joins  the  queue  in  center  2  (with  probability  1-p^). 


(i)  Services  at  centers  1  and  2  are  not  interruptable 

(ii)  Routing  for  type  v  jobs  determined  by  binary  valued  variable  \ 

(iii)  Type  1  jobs  have  non-preemptive  priority  over  type  2  jobs 


Figure  10.1.  Closed  network  of  queues  with  two  job  types 


■-  ? '  </*  V 
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After  completion  of  service  at  center  2,  jobs  join  the  queue  at  center  1. 
At  both  service  centers,  type  1  jobs  have  nonpreemptive  priority  over  type 
2  jobs.  Jobs  of  the  same  type  at  either  of  the  centers  receive  service 
in  order  of  their  arrival  at  the  center.  We  assume  that  all  service  times 
are  mutually  independent;  jobs  of  type  v  at  center  i  receive  service  which 
is  exponentially  distributed  with  parameter  .  The  limiting  response 
time  for  type  v  jobs  that  we  consider  in  this  model  is  the  time 
measured  from  when  upon  completion  of  service  at  center  2,  a  type  v  job 
enters  the  queue  at  center  1,  until  the  next  such  entrance  by  the  job  into 
the  queue  at  center  1. 

In  this  model,  there  are  two  job  classes,  class  1  jobs  at  center  1 
and  class  2  jobs  at  center  2.  Each  center  sees  both  job  types,  but  only 
one  job  class.  The  irreducible  Markov  routing  matrices  are  of  the 

form 


Since  type  1  jobs  have  priority  over  type  2  jobs  at  both  centers,  the 
(type,  class)  pairs  ordered  by  decreasing  priority  are  J^(i)*(l,i)  and 
J2(i)*(2,i),  1*1,2.  For  this  model,  it  is  sufficient  to  take  as  the 
component  S^t)  in  the  vector  Z(t)  the  type  of  job  in  service  at  center  i 
at  time  t,  rather  than  the  (type,  class)  pair.  Then  we  can  define  the 
vector  Z(t)  as 

2<t>  ■  (cf)(t>,C<1>(t),S1(t),C<2>(.),C<1>(t),S2(t))  . 
where,  for  1*1,2  and  v*l,2, 
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C^V\t)  -  number  of  type  v  jobs  In  queue  at  center  1  at  time  t  , 

and 

S1(t)  -  type  of  job  in  service  at  center  1  at  time  t 
0  if  center  1  Is  idle  at  time  t  . 

Letting  Nv(t)  (v-1,2)  denote  the  position  from  the  top  of  the  type  v  marked 
job  in  the  linear  job  stack,  for  t^O  the  state  vector  for  this  model  is 

X(t)  -  (Z(t),N1(t),N2(t))  . 

Letting  L(t)  denote  the  last  state  visited  by  the  Markov  chain  X«{X(t) :t£0} 
before  jumping  to  X(t) ,  the  vector  V(t)  is 

V(t)  -  (L(t) ,X(t) )  . 

For  N-2  jobs,  the  state  space  E  of  the  process  {X(t):t£0}  has  six 
states  and  is 


E  -  {(0,0, 0,1, 0,1, 2,1),  (0,0, 0,0, 1,2, 1,2),  (0,0, 1,0, 0,2, 1,2)} 

u  {(0,0,2, 0,0, 1,2,1),  (1,0, 1,0, 0,0, 2,1),  (0, 1,2, 0,0, 0,1, 2, ) }  . 

The  subsets  A^  and  A^  of  E  defining  the  start  of  response  times  for 
the  type  1  marked  job  are 

-  {(0,0, 0,1,0, 1,2,1),  (0,0, 2,0, 0,1, 2,1)} 


A^1J  -  {(0,0, 1,0, 0,2, 1,2),  (0,1, 2,0, 0,0, 1,2)}  . 
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(2}  (2^ 

Similarly,  the  subsets  A^  and  A^of  E  defining  the  start  of  response 
times  for  the  type  2  marked  job  are 

aJ2)  -  t(0, 0,0,0, 1,2, 1,2),  (0,0, 1,0,0, 2, 1,2)} 

and 

A<2>  -  {(0,0, 2, 0,0, 1,2,1),  (1,0,1, 0,0,0, 2,1)}  . 

Since  and  R^  are  response  times,  b|V^*A^V^  and  B^^A^,  V"l,2; 
see  Figure  10.2.  It  is  easy  to  check  that  the  state  space  F  of  the  process 
{V(t):t20}  has  nine  states.  The  subsets  S^  of  F  defining  the  starts  of 
response  times  for  the  type  \)  marked  job  are 

S(1)  -  {(0,0, 0,1,0, 1,2, 1,0,0, 1,0,0, 2, 1,2)} 

U  {(0,0, 2, 0,0, 1,2, 1,0, 1,2, 0,0, 0,1, 2)} 

and 

S(2)  -  {(0,0, 0,0, 1,2, 1,2, 0,0, 2,0, 0,1, 2,1)} 

u  {(0,0, 1,0, 0,2, 1,2, 1,0, 1,0, 0,0, 2,1)}  , 

respectively;  see  Figure  10.3.  Here  ve  use  the  enumeration  of  the  six 
states  of  £  given  in  Figure  10.2.  Thus,  e.g.,  (1,3)  denotes  the  state 
(0, 0,0, 1,0, 1,2, 1,0,0, 1,0,0, 2, 1,2)«F. 

Simulation  results  for  this  model  for  p^«p^*p*0.75,  x£2^-X^2^«Xj*l 
and  X2^*X22^*X2"0.5,  with  N»2,  appear  in  Tables  10.1-10.4.  With  these 
parameter  values,  there  is  one  type  1  job  and  one  type  2  job.  The  routing 
and  service  requirements  of  the  two  job  types  are  the  same;  the  two  jobs 
differ  only  with  respect  to  the  nonpreemptlve  priority  given  (at  each 
center)  to  the  type  1  job.  The  simulation  used  the  congruentlal  uniform 


Figure  10.2.  State  transitions  in  Markov  chain  X  and  subsets  of  E 
for  response  times  R(1*  and  R<2* 
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r 


Figure  10.3.  State  transitions  in  Markov  chain  V  and  subsets  of  F 
for  response  times  R*1*  and  R{2) 


random  number  generator  described  by  Lewis,  Goodman,  and  Miller  (1969), 
with  exponential  service  times  obtained  by  logarithmic  transformation  of 
the  uniform  random  numbers.  Independent  streams  of  exponential  random 
numbers  (obtained  from  different  seeds)  were  used  to  generate  individual 
exponential  holding  time  sequences. 

For  the  simulation  results  of  Tables  10.1-10.4,  the  return  state 
defining  0-cycles  of  the  response  time  for  the  type  1  job  is  the  state 
(0, 0,2, 0,0, 1,2, 1,0, 1,2, 0,0, 0,1, 2) .  This  corresponds  to  a  response  time 
for  the  type  1  (marked)  job  starting  when  the  type  2  (marked)  job  is  in 
service  at  center  1.  Table  10.1  summarizes  results  of  the  simulation  and 
reports  point  estimates  and  90  percent  confidence  intervals  for  the 
quantities  E{R^},  E{R^}  and  E{R^  }-E{R^ }  over  a  range  of  number  of 
cycles  of  the  type  1  marked  job.  Theoretical  values  for  these  quantities 
are  shown  in  parentheses.  Thus,  for  example,  100  cycles  of  the  type  1 
marked  job  were  observed  in  the  simulated  time  interval  (0,903.00)  and 
there  were  a  total  of  446  transitions  in  the  continuous  time  Markov  chain 
{L(t):t2t0}.  A  total  of  130  response  times  for  the  type  1  (marked)  job 
were  observed  along  with  56  response  times  for  the  type  2  (marked)  job. 
For  the  quantity  E{R^}-7,  the  point  estimate  6.946  was  obtained,  and 
the  90  percent  confidence  interval  had  half-length  0.6334.  Note  that  for 
E{R^)}  and  E{R^},  all  of  the  confidence  intervals  surround  the 
theoretical  values.  In  the  case  of  E{R^ }-E{R^ } ,  the  confidence 
intervals  based  on  Equation  (10.2.8)  also  surround  the  theoretical  value. 
Table  10.2  gives  results  obtained  for  P{R^£x),  with  x*4,  8,  12,  16  and 
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In  Table  10.3  we  give,  for  the  several  values  of  x ,  point  and  Interval 
estimates  for  P{R^£x}-P{R^£x} ,  based  on  the  use  of  super-cycles  and 
Equation  (10.2.12).  Thus,  for  x* 4,  100  cycles  based  on  response  times  for 
the  type  1  job  resulted  In  37  super-cycles  defined  by  response  times  for 
the  type  2  job  greater  than  x.  Note  that  the  number  of  cycles  for  the 
type  1  marked  job  has  been  fixed,  and  for  each  x  the  estimates  for 
P{R^£x}-P{R^£x)  computed  from  the  resulting  random  number  of 
super-cycles . 

Table  10.4  contains  estimates  of  the  quantities  P{R'  1  £x}  obtained 
from  the  standard  regenerative  method  applied  to  these  super-cycles.  An 
overall  observation  from  Tables  10.2  and  10.4  Is  that  the  lengths  of 
confidence  intervals  obtained  for  P{R^£x}  and  P{R^£x}  are  roughly 
comparable. 
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TABLE  10.2 

Percentiles  of  Type  1  Response  Times  in  Closed  Network  of  Queues 
With  Two  Job  Types.  N^-l,  N2“l,  p“0.75,  Xj_«l,  X2*0.5. 
Return  State  is  (0, 0, 2,0, 0,1, 2, 1,0, 1,2, 0,0, 0,1, 2) . 
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TABLE  10.3 

Difference  of  Percentiles  of  Response  Times  In  Closed  Network  of  Queues 
With  Two  Job  Types.  N^*l,  N2*l,  p-0.75,  A^-l,  \2"0.5. 

Return  State  is  (0,0,2, 0,0, 1,2, 1,0,1, 2, 0,0,0, 1,2) . 


P{R(1)s;4}-P{R(2)s4} 
No.  of  super-cycles 

P{R(1)s8}-P{R(2^8} 
No.  of  super-cycles 


P{R(1)S12}-P{R(2)S12} 
No.  of  super-cycles 

P{R(1)S16}-P{R(2)S16} 
No.  of  super-cycles 


P{R(1)S20}-P{R(2)S20} 
No.  of  super-cycles 


No.  of  Cycles  For  Type 

1  Marked  Job 

100 

200 

400 

800 

1000 

0.1254 

0.1295 

0.1342 

0.1332 

0.1384 

±0.0784 

±0.0627 

±0.0417 

±0.0275 

±0.0192 

37 

79 

181 

347 

438 

0.3988 

0.3409 

0.3111 

0.3036 

0.3125 

±0.1227 

±0.0817 

±0.0525 

±0.0373 

±0.0201 

29 

61 

134 

254 

319 

0.3915 

0.3543 

0.3331 

0.3259 

0.3280 

±0.1424 

±0.1112 

±0.0677 

±0.0496 

±0.0140 

22 

46 

93 

180 

224 

0.2850 

0.2970 

0.2693 

0.2693 

0.2636 

±0.1288 

±0.1092 

±0.0606 

±0.0489 

±0.0088 

15 

32 

65 

129 

153 

0.2422 

0.2470 

0.2119 

0.2078 

0.2009 

±0.1081 

±0.0825 

±0.0530 

±0.0415 

±0.0043 

11 

24 

43 

84 

104 

170 


TABLE  10.4 

Percentiles  of  Type  2  Response  Times  In  Closed  Networks  of  Queues 
With  Two  Job  Types.  Nj»l,  N2*l,  p“0.75,  Xj»l,  X2“0.5. 

Return  State  is  (0,0, 2, 0,0, 1,2, 1,0,1, 2, 0,0, 0,1, 2) . 
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11.0.  IMPLEMENTATION  CONSIDERATIONS 

In  order  to  carry  out  a  passage  time  simulation  of  a  network  of  queues, 
we  must  be  able  to  generate  sample  paths  or  realizations  of  the  stochastic 
system.  A  necessary  part  of  any  such  generation  procedure  is  an  algorithm 
(or  algorithms)  for  random  number  generation,  i.e.,  for  the  generation  of 
numbers  that  can  be  treated  as  instances  (samples)  of  random  variables. 

In  this  section  we  consider  aspects  of  random  number  generation  pertinent 
to  the  implementation  of  passage  time  simulations  according  to  the  methods 
of  the  previous  sections. 

11.1.  Random  Number  Generators 

Our  discussion  follows  Learmonth  and  Lewis  (1973a) .  By  a  "random 
number  generator"  (or  "pseudo-random  number  generator")  we  mean  an 
algorithm  which  produces  sequences  of  numbers  that  follow  a  specified 
probability  distribution  and  possess  the  appearance  of  randomness.  The 
use  of  "sequence  of  numbers"  means  that  the  algorithm  is  to  produce  many 
random  numbers  in  a  serial  fashion.  Even  though  a  particular  user  may 
need  only  relatively  few  of  the  numbers,  we  generally  require  that  the 
algorithm  be  capable  of  producing  many  numbers.  "Probability  distribution" 
implies  that  we  can  associate  a  probability  statement  with  the  occurrence 
of  each  number  produced  by  the  algorithm.  We  usually  take  the  probability 
distribution  to  the  uniform  distribution  on  the  interval  [0,1].  If  a 
source  of  [0,1]  uniform  random  numbers  is  available,  then  in  principle  it 
is  possible  to  transform  these  uniform  random  numbers  by  means  of  the 
Inverse  probability  integral  into  random  numbers  having  any  desired 
distribution.  For  reasons  of  computational  efficiency,  however,  a  large 
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amount  of  effort  has  gone  Into  the  development  of  methods  for  direct 
generation  of  random  numbers  having  nonuniform  distributions;  see  Ahrens 
and  Dieter  (1973a)  for  a  comprehensive  discussion.  With  respect  to 
"appearance  of  randomness,"  it  may  be  somewhat  surprising  that  the  actual 
implementation  of  most  commonly  used  algorithms  for  uniform  random  number 
generation  is  as  a  (deterministic)  recurrence  relation  in  which  each 
succeeding  number  is  a  function  of  the  preceding  number.  Thus,  although 
true  randomness  requires  independence  of  successive  numbers,  the  algorithm 
generates  a  deterministic  dependent  sequence.  When  parameters  of  the 
recurrence  relation  are  chosen  carefully,  such  algorithms  for  uniform 
random  number  generation  do  yield  sequences  which  (statistically)  appear 
to  be  random.  This  appearance  of  randomness  is  the  origin  of  the  term 
"pseudo-random  numbers . " 

Since  the  results  of  a  simulation  depend  critically  on  an  acceptable 
appearance  of  randomness,  it  is  important  that  a  proposed  uniform  random 
number  generator  be  subjected  to  thorough  statistical  testing.  Although 
che  simulation  practitioner  need  not  necessarily  be  concerned  with  the 
details  of  the  rather  specialized  techniques  for  statistical  testing  of 
random  number  generators,  he  should  be  convinced  prior  to  use  that  an 
available  uniform  random  number  generator  has  been  successfully  tested. 

See  Fishman  (1978),  Ch.  8  for  a  discussion  of  statistical  tests  for  uniform 
random  number  generators. 

The  most  widely  used  (uniform)  random  number  generators  are  of  a  class 
•  iwM  a#  linear  congruential  generators.  Such  generators  employ  a 
•  .r 'Mice  relation  of  the  form 
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XQ  -  bXQ_^  +  c  (mod  m)  (11.1.1) 

In  Equation  (11.1.1)  all  quantities  are  nonnegative  Integers.  This 

equation,  read  "X  equals  bX  ,+c  modulo  m,"  says  that  X  is  the  remainder 
n  n-i  n 

when  bXn_^+c  is  divided  by  m.  The  quantity  b  is  called  the  multiplier. 

m  is  the  modulus,  and  c  is  the  increment.  Given  a  startina  value  v° 

and  values  b£0,  c2Q ,  and  m  such  that  m>XQ,  m>b  and  m>c,  a  sequence  of 

integers  X^.X^,..*  is  generated  by  successive  application  of 

Equation  (11.1.1).  Uniform  random  numbers  U  on  the  interval  [0,1]  are 

n 

obtained  by  dividing  by  m,  l.e.,  for  n~l,2,..., 

U  -  X  /m  (11.1.2) 

n  n 

The  recurrence  relation  of  Equation  (11.1.1)  is  sometimes  called  a  "mixed 
linear  congruential  generator,”  the  term  "mixed"  coming  from  the  fact  that 
it  involves  a  multiplication  by  a  constant  b  along  with  an  addition  of  a 
constant  c.  Many  random  number  generators  are  "multiplicative"  or  "pure 
congruential"  in  that  c-0,  giving 

X  -  bX  .  (mod  m)  .  (11.1.3) 

n  n-i 

The  initial  or  starting  value  XQ  is  often  called  the  seed  of  the  random 
number  generator. 

Although  it  may  appear  that  Equation  (11.1.1)  produces  m  distinct 
numbers,  this  is  not  the  case  unless  b  and  m  are  chosen  properly.  It  is 
characteristic  of  generators  of  this  type  that  there  is  ultimately  a  cycle 
of  numbers  which  is  repeated  indefinitely;  this  repeating  cycle  of  numbers 
is  called  the  period  of  the  generator.  It  is  clear  that  a  congruential 
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sequence  used  as  a  source  of  random  numbers  should  have  a  long  period, 
and  since  the  period  can  never  be  greater  than  m,  the  value  of  m  should 
be  rather  large. 

Mathematical  results  based  on  number- theoretic  considerations  are 
available  for  characterizing  the  values  of  b  and  c  which  result  in  the 
maximum  period  length  m;  see  Knuth  (1969),  Ch.  3.  For  the  special  case  of 
multiplicative  congruential  generators  (c*0) ,  the  basic  result  concerning 
maximum  period  length  says  that  the  maximum  period  length  (m)  is  not 
achievable.  It  is,  however,  still  possible  to  obtain  multiplicative 
congruential  generators  with  quite  long  periods.  Results  characterizing 
the  maximum  period  for  this  multiplicative  case  are  available,  but  the 
number-theoretic  considerations  are  Involved. 

If  the  modulus  m  in  a  multiplicative  congruential  generator  is  prime, 
(i.e.,  has  no  divisors  other  than  1  and  itself)  a  period  of  length  m-1  is 
achievable.  Such  a  period  length,  of  course,  is  just  one  less  than  the 
maximum  possible  length.  If,  in  addition,  we  choose  the  multiplier  b  so 
as  to  satisfy  an  appropriate  (sufficient)  number-theoretic  condition  with 
respect  to  (prime)  m,  then  for  any  starting  value  Xg<m,  the  maximum  period 
length  m-1  is  achieved.  The  determination  of  values  for  multipliers  b 
satisfying  the  number-theoretic  condition  for  maximum  period  length  in  a 
multiplicative  congruential  generator  in  general  Involves  lengthy 
calculations.  Further  details  are  in  Knuth  (1969),  Ch.  3. 
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Zn  any  particular  digital  computer  system,  only  a  finite  number  of 

positive  integers  are  representable,  the  limitation  being  the  word  size 

of  the  system.  We  now  state  a  particular  (multiplicative  congruentlal) 

uniform  random  number  generator  which  utilizes  the  full  word  size  of  IBM 

System/360  (370)  computer  systems.  (This  is  the  uniform  random  generator 

used  to  obtain  the  numerical  results  in  previous  sections.)  In  the 

System/360,  the  word  size  is  32  bits  with  1  bit  reserved  for  algebraic  sign; 

31 

an  obvious  choice  for  m  is  thus  2  .  A  multiplicative  congruentlal 

k 

generator  with  m“2  (for  some  positive  Integer  k)  can  have  a  maximum  period 

31 

length  of  m/4.  Thus  for  System/ 360  computer  systems  with  m-2  ,  the 

29 

maximum  period  is  2  ,  and  the  period  length  may  also  depend  on  the 

starting  value.  It  happens  (fortuitously)  that  the  largest  prime  less 

31  31  31 

than  or  equal  to  2  is  2  -1.  Hence,  by  choosing  m«2  -1,  it  is  possible 

to  implement  on  System/ 360  computer  systems  uniform  random  number 

31 

generators  having  a  maximum  period  length  of  2  -2.  Note  that  the 

number-theoretic  conditions  ensuring  a  maximum  period  length  do  not 
necessarily  guarantee  good  statistical  properties  for  the  generator, 
although  the  choice  of  the  particular  multiplier  73  does  satisfy  some 
known  conditions  regarding  statistical  properties  of  the  resulting 
sequence. 

System/ 360  Generator 

Let  Xq>0.  Then  for  n^l, 

X„  -  75  X  .  (mod  231-1) 
n  ti“i 

-  16807  X  .  (mod  231-1)  (11.1.4) 

n~i 


f 
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and 

U  -  X  / (231-1)  (11.1.5) 

n  n 

Tha  uniform  random  numbar  ganarator  of  Equation  (11.1.5)  has  bean  tested 

extensively,  and  the  results  of  the  statistical  tests  Indicate  that  it  is 

very  satisfactory;  see  Lewis,  Goodman  and  Miller  (1969)  and  Learmonth  and 

31 

Lewis  (1973b) .  Other  multipliers  for  generators  with  modulus  m«2  -1  are 

in  use.  Results  of  pertinent  statistical  tests  are  given  by  Hoaglin 
(1976). 

11.2.  Nonuniform  Random  Numbers 

The  problem  of  generating  random  numbers  from  a  specified  (nonuniform) 
distribution  is  in  principle  solved  by  having  a  source  of  uniform  random 
numbers  and  transforming  these  random  numbers  by  means  of  the  inverse 
probability  integral.  Because  it  is  not  always  possible  to  compute  or  to 
compute  efficiently  the  inverse  of  a  given  distribution  function,  a  great 
deal  of  effort  has  gone  into  the  development  of  methods  for  direct 
generation  of  nonuniform  random  numbers;  see  Ahrens  and  Dieter  (1973b) 
for  a  comprehensive  discussion.  Desirable  properties  of  such  direct 
methods  are  that  they  be  exact,  very  fast,  and  economical  of  computer 
storage.  The  property  of  exactness  is  that  any  deviation  from  the 
specified  distribution  results  from  computer  round-off  error  rather  than 
a  defect  in  the  method  itself.  Comparisons  are  hard  to  make  among 
particular  methods,  partly  because  of  machine  dependencies.  It  is, 
however,  almost  always  true  that  with  very  little  cost  in  complexity,  it 
is  possible  to  improve  on  the  Inverse  probability  Integral  transformation 

i 
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by  an  ordar  of  magnitude  in  execution  time;  the  fastest  available 
algorithms  for  nonuniform  random  number  generation,  however,  are  often 
quite  complex. 

The  basis  for  the  generation  of  nonuniform  random  numbers  by 
transformation  of  a  uniform  random  number  is  the  following  statement.  If 
U  is  uniformly  distributed  on  [0,1]  and  if  F(x)  is  any  distribution 
function,  then  the  random  variable 

X  -  F_1(U) 

has  distribution  F(x);  for  OSuil,  the  inverse  function  F_1(u)  is  defined 
by 

F-1(u)  *  inf {z:F(z)Su}  . 

It  follows  that  to  generate  samples  of  a  random  variable  X  having 
distribution  F(x)  from  a  random  number  U  (uniformly  distributed  on  [0,1]), 
we  must  be  able  to  solve  the  equation 

F(x)  »  u  . 

Then  given  a  uniform  random  number  U,  we  return 

X  •  inf {z:F(z)*U}  . 

Figure  11.1  illustrates  the  inverse  transformation  method.  Note  that 
this  technique  applies  to  discrete  as  well  as  continuous  random  variables. 


The  Inverse  transformation  method  provides  a  straightforward  means  of 
generating  samples  of  an  exponential  random  variable,  as  may  be  needed 
for  a  passage  time  simulation  of  a  network  of  queues.  We  can  obtain  an 
exponential  (rate  parameter  X)  random  number  X  by  generating  U,  a  uniform 
random  number  on  [0,1],  and  transforming  it  according  to 

X  -  -(£n  U)/X  . 

This  transformation  is  obtained  by  solving  the  equation  U*»F (X)  for  X, 
yielding  X*{-£n(l-U)}/X,  and  observing  that  1-U  is  uniformly  distributed 
on  the  interval  [0,1]  when  U  is.  Note,  however,  that  although  this 
logarithmic  transformation  is  easy  to  implement,  it  is  a  relatively  slow 
method  for  obtaining  exponential  random  numbers.  The  fastest  methods 
available  at  present  for  exponential  random  numbers  use  so-called 
"decomposition"  methods  using  the  ideas  of  Marsaglia,  MacLaren  and  Bray 
(1964) .  The  basis  for  the  method  is  division  of  the  random  variable  into 
several  populations,  from  most  of  which  samples  can  be  obtained  easily. 
Using  geometric  considerations,  the  density  function  of  the  exponential 
random  variable  is  decomposed  into  a  large  number  of  rectangular  regions, 
wedge-shaped  regions,  and  a  tail.  The  Navel  Postgraduate  School  random 
number  generator  package  LLRANDOM  (Learmonth  and  Lewis  (1973a))  contains 
an  IBM  System/360  Basic  Assembler  Language  implementation  of  a 
decomposition  method  for  generation  of  exponential  random  numbers. 

In  implementation  of  a  passage  time  simulation,  it  may  be  necessary 
to  generate  samples  of  a  random  variable  X  having  a  Cox-phase  (exponential 
stage)  representation.  It  is  easy  to  do  so  if  generators  of  uniform  and 
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exponential  random  numbers  are  available.  Using  the  notation  o£ 

Section  3.1,  suppose  that  the  distribution  of  X  has  n  exponential  stages. 
For  j-l,2,...,n, 


X  -  xx  +  x2  +  . . .  +  Xj 

with  probability  (l-bj)aj,  where  the  random  variables  X^ ,Xg » . . . ,X^  are 
mutually  independent  and  X^  is  exponentially  distributed  with  rate 
parameter  X^.  To  obtain  random  samples  of  X,  we  use  a  two  step  procedure. 
First,  we  generate  U,  a  uniform  random  number  on  [0,1],  and  set  k  equal  to 
1  if  CKUSd-b^Jaj^  and  equal  to 

min{  j  >1 :  (1-b ^  _1)  a^  (1-b ^  )  a^  } 

otherwise.  Then,  using  the  inverse  transformation  method  for  exponential 
random  numbers,  we  generate  k  (mutually  Independent)  uniform  random 
numbers  Ui,U2**<*,Uk  and  return 

k 

X  -  S  (-in  U.)/X  . 

j-1  J  J 

The  inverse  transformation  method  for  generating  random  numbers  having 
a  specified  discrete  distribution  provides  a  means  of  routing  jobs  through 
a  network  of  queues.  In  the  simplest  case  (e.g.,  in  the  network  of 
Figure  S.l),  it  is  necessary  to  generate  samples  of  a  Bernoulli  random 
variable  W  for  which  (with  0<p<l) 

P{W-1}  -  p 

and 

P{W-0}  -  1-p  . 
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To  do  so,  we  generate  a  random  number  U  (uniformly  distributed  on  [0,1]), 
and  return  W-l  if  U£p,  and  return  W*Q  otherwise.  By  this  procedure,  we 
are  in  effect  partitioning  the  interval  [0,1]  and  determining  the  value 
of  X  by  the  portion  of  the  interval  in  which  the  generated  uniform  random 
number  lies.  The  generalization  to  handle  more  complex  routing  from  a 
service  center  is  straightforward. 

11.3.  Single  and  Multiple  Streams  of  Random  Numbers 

It  is  typically  the  case  when  simulating  a  network  of  queues  that  we 
have  need  for  several  streams  of  random  numbers  (e.g.,  random  service 
times  for  each  of  several  servers,  routing  through  the  network,  etc.). 

Since  most  algorithms  for  nonuniform  random  number  generation  require  the 
generation  of  one  of  more  uniform  random  numbers,  the  question  arises  as 
to  whether  single  or  multiple  streams  of  uniform  random  numbers  should  be 
used.  When  a  single  stream  of  uniform  random  numbers  is  used,  the  course 
of  the  simulation  determines,  usually  in  a  complex  manner,  the  role  of 
individual  uniform  random  numbers;  thus,  e.g.,  for  the  cyclic  queues  of 
Figure  2.1,  a  random  subsequence  of  the  generated  uniform  random  numbers 
can  be  transformed  to  give  the  service  times  at  one  of  the  service  centers, 
with  the  remaining  random  numbers  used  to  generate  the  service  times  at 
the  other  center.  Alternatively,  if  appropriate  seeds  are  available,  we 
can  use  nonoverlapping  portions  of  the  uniform  random  number  sequence  to 
generate  the  service  times  at  the  individual  centers.  The  concern  is  that 
when  a  single  stream  of  uniform  random  numbers  is  used,  we  are  in  effect 
assuming  that  particular  random  subsequences  of  the  original  uniform  random 
number  sequence  has  an  acceptable  appearance  of  randomness,  and  this  may 
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not  be  the  case.  There  ere  examples  of  simulations  where  the  use  of  a 
single  stream  of  uniform  random  numbers  has  led  to  rather  bizarre  results. 
Although  this  aspect  of  random  number  generation  is  not  well -understood, 
in  many  cases  it  is  probably  good  practice  to  use  separate  streams.  Some 
additional  bookkeeping  is  of  course  necessary  to  handle  the  separate 
streams,  and  judgement  is  required  as  to  what  extent  multiple  streams 
should  be  used  when  simulating  a  complex  stochastic  system.  For  a  network 
of  queues,  it  is  probably  advisable  to  use  separate  random  number  streams 
for  the  interarrival  times,  service  times,  and  routing  of  jobs  from  the 
individual  service  centers. 

In  Table  11.1  we  give  values  of  seeds  which  can  be  used  to  generate 
Independent  streams  of  uniform  random  numbers  from  the  System/360  uniform 
random  number  generator  of  Equations  (11.1.4)  and  (11.1.5).  These  seeds 
are  values  of  X  which  are  100000  apart  in  the  sequence  of 

Cl 

Equation  (11.1.4),  i.e.,  if  Xq-377003613,  then  x100000"’648473574’ 

X2QQQQQ» 1396717869,  etc.  It  is  necessary  when  using  multiple  streams  of 
random  numbers  to  keep  in  mind  approximately  how  many  random  numbers  are 
needed;  undesired  dependence  among  random  numbers  may  result  if  portions 
of  the  original  sequence  overlap  inadvertently. 

11.4.  Generation  of  State  Vector  Frocesses 

When  carrying  out  a  simulation  of  a  network  of  queues  (or  any 
stochastic  system) ,  we  observe  the  behavior  of  the  system  as  it  evolves 
in  time.  Implicit  in  any  implementation  of  the  simulation  is  the 
definition  of  an  appropriate  system  state  vector.  This  "state  of  the 
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system  at  time  t"  constitutes  a  stochastic  process,  and  to  obtain  estimates 
of  quantities  of  interest,  we  must  somehow  generate  realizations  or  sample 
paths  of  this  state  vector  process.  For  complex  networks  of  queues,  it 
is  often  convenient  to  generate  the  process  (e.g.,  using  an  event 
scheduling  approach)  by  means  of  timing  routines  applicable  to  a  general 
discrete  event  simulation,  as  typically  provided  by  a  high  level  simulation 
programming  language.  However,  when  there  is  a  characterization  of  the 
state  vector  process  as  a  familiar  stochastic  process,  it  may  be  possible 
to  generate  the  process  directly  and  more  efficiently  (with  respect  to 
speed)  than  by  using  all  of  the  apparatus  for  timing  which  is  necessary 
for  a  general  discrete  event  simulation.  This  is  relevant  to  the  passage 
time  simulations  discussed  here  in  that  the  state  vector  process  {X(t):t^0} 
of  Section  3.2  constitutes  a  finite  state  continuous  time  Markov  chain. 

We  consider  generation  of  such  Markov  chains  next. 

Let  Y«{Y(t):t20}  be  a  continuous  time  Markov  chain  having  finite  state 
space  E,  and  let 


a  - 

be  its  matrix  of  Infinitesimal  transition  parameters;  thus  for  i,jeE, 
q1jkO  for  ij*j,  and  for  all  i, 


Denote  by  {x^inSO}  the  jump  times  of  the  process,  and  for  n-0,1,...,  set 

Y  -  Y(x  )  . 
n  n 


\ 

•  i 
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Generation  of  the  continuous  time  Markov  chain  %  can  be  based  on  the 
following  characterization  (see,  e.g.,  ^inlar  (1975a),  p.  247).  For  any 
jeE,  ueR+,  and  n-0,1,..., 

P*Yn+l"^’  Ttt+l”Tn>u^Y0 . Yn;  T0 . Tn^  "  rije  ***  (11.4.1) 

if  XQ«*i.  Here,  q^- q^,  r^-q^^  for  ^  •  and  rii“0,  Thu8»  given  a 
jump  to  state  1,  the  process  remains  in  state  i  for  an  exponentially 
distributed  (rate  parameter  q^)  amount  of  time,  and  then  jumps  to  state  j 
with  independent  probability  r^.  It  follows  that  generation  of  a  sample 
path  for  the  Markov  chain  Y,  (i.e.,  the  sequence  of  jump  times  and 
successive  states)  can  be  accomplished  by  successive  generation  of  a  pair 
of  independent  random  numbers.  This  pair  consists  of  an  exponential  random 
number  and  a  sample  from  a  discrete  distribution  specified  by  the  jump 
probabilities.  Note  that  each  element  of  such  a  pair  can  be  generated  by 
the  inverse  transformation  method. 

This  discussion  of  the  generation  of  continuous  time  Markov  chains 
presupposes  that  the  elements  of  the  ^.-matrix  (and  hence  the  jump 
probabilities)  are  available  explicitly.  For  the  class  of  networks  of 

mJ 

queues  discussed  here,  enumeration  of  the  state  space  and  explicit 
calculation  of  the  Infinitesimal  transition  parameters  is  in  general 
somewhat  tedious.  It  is  important  to  observe,  however,  that  complete 
knowledge  of  the  ^-matrix  (e.g.,  for  the  Markov  chain  defined  by 
Equation  (3.2.3))  is  contained  in  the  (given)  routing  matrix  1?  and  the 
parameters  of  the  (exponential  or  Cox-phase)  service  times.  Consequently, 
based  on  Equation  (11.4.1),  it  is  often  possible  to  construct  a  more 
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efficient  algorithm  for  direct  simulation  of  a  given  network  of  queues 
than  that  which  results  from  a  general  discrete  event  simulation.  For 
the  DL/I  component  model  of  Section  5.2,  Appendix  4  gives  such  an  algorithm. 


/ 


'Tf  **■**-. 
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TABLE  11.1 

Seads  for  System/360  Uniform  Random  Number  Generator. 

Values  of  X  100,000  apart  in  X  -75  X  -1  (mod  231-1) . 
n  q  n 

(to  be  read  across) 


377003613 

648473574 

1396717879 

2027350275 

1356162430 

1752629996 

745806097 

201331468 

1393552473 

1966641861 

711072531 

769795447 

1074543187 

1933483444 

625102656 

1116874679 

1442211901 

989455196 

1996695068 

1850124212 

1267310126 

1741371275 

886499692 

1014119573 

933913228 

2082204497 

920168983 

1079618777 

1888797415 

1002901030 

1582733583 

254293472 

1095895189 

219529399 

1706847402 

1951007719 

1169002398 

1482199345 

1976077334 

775245191 

1976418161 

35067978 

400884188 

1895732964 

1904749580 

1301700180 

63685808 

936615625 

110322717 

1029730003 

251900732 

725094089 

828842333 

1471230052 

1703522097 

1356420548 

1670372925 

437765009 

39279049 

2123613511 

150006407 

1633650593 

751601611 

1410990605 

1262214427 

645360044 

1504645702 

1063375004 

941885586 

1753135176 

253642018 

1701685042 

1448665492 

1034856864 

428280431 

259758456 

600732272 

704726097 

398944698 

114386769 

288727775 

1499601820 

2136214308 

1197972807 

1888007825 

686553263 

747119178 

154337000 

136758808 

9182540 

303111010 

154232008 

921093990 

1684263351 

1166344707 

1167753617 

1374693082 

1812641667 

502455872 

857532898 
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APPENDIX  1.  CONVERGENCE  OF  PASSAGE  TIMES 

Label  the  jobs  from  1  to  N  end  for  tiO  set 

Y(t)  -  (Z(t)lS1(t),S2(t) . N^t))  , 

where  N*(t)  denotes  the  position  In  the  job  stack  et  time  t  of  the  job 
labelled  1.  The  vector  Z(t)  is  the  sane  as  in  Section  3.2.  Also  define 
the  marginal  processes 

yV)  -  (ZCO^Ct)) 

for  i»l,2 . N.  All  of  the  processes  (Y(t) :t£0}  and  {Y*(t):t^O}  are 

irreducible,  positive  recurrent  continuous  time  Markov  chains  defined  on 
a  common  underlying  probability  triple,  (Q^P) ,  say.  Observe  that  if  the 
marked  job  of  Section  3  is  the  job  labelled  i,  then  the  process  {Y*(t):tkO} 
coincides  with  the  process  {X(t):tiO}  defined  by  Equation  (3.2.3),  except 
possibly  for  the  Initial  condition  at  t“0.  Define  for  each  job  two 
sequences  of  times,  the  starts  and  terminations  of  the  successive  passage 
times  for  the  job.  For  the  job  labelled  i,  denote  these  times  by  {S*:j20} 

J 

and  {T^:j2l}.  The  definition  of  these  times  in  terms  of  the  process 
{Y*(t):t20}  is  completely  analogous  to  what  was  done  in  Section  4.2  in 
terms  of  the  process  {X(t):t20}.  Then  the  jth  passage  time  for  the  job 
labelled  i  is  Also  define  Markov  chains  (XjijaO)  for 

each  job  in  which  X*  denotes  the  state  of  the  Markov  chain  (Y*(t) stiO) 

J 

when  the  (j+l)8t  passage  time  starts  for  job  i:  X^-Y^S*).  At  this  point 
we  have  N  Markov  renewal  processes,  { (Xj ,S^) :j^0},  all  defined  on  (fl,#?). 
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Next  we  introduce  e  new  sequence  of  passage  times,  {Pj:j£l},  also 
defined  on  (fl,#P);  this  is  the  sequence  of  passage  times  irrespective  of 
job  identity,  enumerated  in  order  of  start  times.  For  each  j,  Pj  is  a 
random  member  of  the  set  {pjsisisj};  this  means  that  Pj»p£|j ^ ,  where  2(j) 
and  k(j)  are  random  variables. 


The  principal  result  of  this  Appendix  is  to  show  that  all  of  the 
sequences  {P^rj^l}  and  {P* : j£l}  converge  in  distribution  to  a  common  random 


variable  P. 


(Al.l)  PROPOSITION.  For  i-l,2,...,N,  P*->P  as  j-*».  In  addition,  Pj->P 
as  j-*“. 

Proof.  Since  the  N  jobs  are  Identical  with  respect  to  their  service 
requirements  and  branching  probabilities,  the  semi-Harkov  kernels  governing 
the  Markov  renewal  processes  { (X*,S*) : jkO)  all  coincide  with  the  kernel  K 

J  J 

of  Section  4.2.  In  fact,  for  any  particular  job  the  only  difference  from 
the  setup  of  Section  4.2  is  that  (with  possibly  one  exception)  the  job  does 
not  start  a  passage  time  at  t-0.  However,  this  difference  does  not  alter 
limiting  results;  the  job  labelled  i  starts  a  passage  time  with  probability 
one  (since  {Y*(t):t20}  is  positive  recurrent),  and  once  this  occurs,  the 
situation  is  exactly  as  in  Section  4.2.  Note  in  particular  that  S^-H®  a.s. 
for  all  1;  thus,  there  is  always  a  next  passage  time  for  every  job.  This 
being  so,  we  have  Pj«">P  as  j**»,  just  as  was  the  case  in  Section  4.2.  (This 
result  is  to  be  expected  since  the  marked  job  was  selected  arbitrarily.) 
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Next  v«  show  that  Pj*oP.  sines  Pj“>P  for  all  1,  wa  can  usa  Cha 
Skorohod  raprasancatlon  theorem  (saa  Skorohod  (1956)  or  Bllllngslay  (1971)) 
Co  assarc  tha  axlstanca  of  a  probability  spaca  (ftJs^P)  and  random  varlablas 
Pj  (jil,  lsisN)  and  P  defined  on  that  spaca  such  Chat  P*  and  P  have  cha 
same  distributions  as  P^  and  P,  respectively,  and  P^P  a.s.  as  j-*00  for 
all  1.  These  representatives  P*  also  provide  representatives  for  tha  P^ 
which  we  call  P^. 


Putting  aside  the  null  sets  of  ft  on  which  the  above  convergence 
statements  do  not  hold,  we  examine  the  numerical  sequence  {Pj(o>):j£l}  for 
one  of  the  remaining  uieft.  We  use  the  following  criterion  for  convergence 
of  a  numerical  sequence  {x^rjal}:  Xj**tc  as  j-*“  if  and  only  If  for  each 
subsequence  (x^ ,}  there  exists  a  further  subsequence  {x^,,}  that  converges 
to  x;  see  Billingsley  (1968),  p.  15,  for  a  similar  usage  In  weak  convergence 
theory.  Select  a  subsequence  {Pj , (w)  } .  This  subsequence  must  contain  a 
further  subsequence  {P^ ,,  (u) }  that  is  identical  to  a  subsequence  of  one  of 
the  sequences  {P^ , (w) :J21},  say  for  i*ig.  This  follows  from  the  fact  that 
there  are  only  a  finite  number  of  jobs.  But  this  subsequence  IP^ , ) 
converges  a.s.  to  P(w)  since  the  full  sequence  does.  Thus  Pj-*-P  a.s.  and 
therefore  Pj-oP. 
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APPENDIX  2.  PROOF  OF  RATIO  FORMULA 

W«  provide  a  proof  of  ch«  ratio  formula  (Theorem  (4.2.9))  for  E{f(X,P)> 
which  makes  it  possible  to  use  the  regenerative  method  for  estimation  of 
passage  times.  The  proof  given  here  does  not  require  the  key  renewal 
theorem. 

(A2.1)  PROPOSITION.  Assume  that  E{|f (X,P) | }<*.  Then 

E{f (X,P)}«E{Y^(f)}/E{a^},  where  Y^(f)  is  given  by  Equation  (4.2.6). 

Proof .  Assume  fiO,  E{f(X,P)}<*,  and  set  f  ■min(f,c)  for  some  c  such  that 

c 

0<c<“.  Clearly,  Equation  (4.2.8)  also  holds  for  f  ,  i.e., 

c 

£e<VW  ~  '=«-P>  • 

and  since  f  is  bounded , 

lim  E{fc(Xn,Pn+1)}  -  E{f  (X,P) }  .  (A2.2) 

nr1*0 

Next  we  compute  the  Cesaro  average  of  the  sequence  appearing  in 
Equation  (A2.2) .  First  we  write 

\  |  i(m)+l  \ 

£  £c(Xn'Pn+l)  {  /<■«■>  "  E  |  2  VV  j  /(o+1)  ~  E{Y' (m) }/ (nrfl)  ,  (A2.3) 

6i(m)+l 

where  l(m)-max{k:g.  Sm)  end  Y'  (m)-  T,  f  (X.P  ..). 

k  n=5+l  c  n  n+1 

Since  0sf,sc,  we  have  OSY' (m)Sc(84(n)+1-o) .  In  addition,  Wald's  equation 
(see  Chung  (1968),  Theorem  5.5.3)  implies  that 
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E{8Jl(m)+l}  "  *(«1>B{tW+l> 
and 

j  *<»)+!  j 

E|  £  Bk(*c>j  "  ^(gJBUW+U  . 

That*  aquations  plus  tha  alamantary  ranawal  theorem  (Smith  (1958),  p.  246) 


imply  that 


lia  E{81(B)+1-m)/(»fl)  -  0 


j Mm)+1  | 

limE)  £  Y.(f  )(  /(m+1)  -  E{Y.  (f  ) }/E{a. }  . 

l  k-1  16  c  f  1  c  1 


Banca  from  Equation  (A2.3)  wa  have 


11m  E 


ijww  j  /<*»  -  «w>  /'<«!>  • 


(A2.4) 


From  Equations  (A2.2)  and  (A2.4)  wa  concluda  that 


E{f  (X,P)>  -  E{Y.(f  )}/E{a-}  . 

C  1  C  X 


(A2.5) 


Now  wa  lat  e*»  on  both  sidas  of  Equation  (A2.5)  and  usa  tha  assumption  that 
E{f «,?)}<•  to  obtain 


E{f(X,P)}  -  E{Y1(f))/E{d1)  . 


(A2.6) 


For  a  general  f  function,  wa  wrlta  f«f  -f"  and  apply  tha  abova  argument  to 
both  f*  and  f”.  Thua  wa  hava  Equation  (A2.6)  provided  E{|f(X,P) |}<“. 


APPENDIX  3.  ESTIMATION  OF  VARIANCE  CONSTANT'S 

2 

W«  first  consldsr  estimation  of  ths  variance  constant  o  appearing  In 
Equation  (10.2.8)  which  loads  to  a  confidancs  interval  for  E(R^ }-E{R^ } . 
Based  on  n  cycles,  for  1*1,2,  compute  as  an  estimate  of 

aii  *  *  vertex^}  -  2r (1)cov{ok,N^1^}  +  (r^^variN*1*) 


according  to 


8ii  •  ‘u 


»(«,(«  +  «(«,*,«> 

"n  *12  +  Un  }  *22  * 


where 


n 

*11  •  2(#J-“n)2  > 

•g>  -  (n-l)"1 


5(1)  x 
\  )  * 


.g>  -  (n-l)”1  23 (M^1>-K1J1> ) 2  , 

i 

a  -  n”1  £a.,  H*1*  -  n*1  23n.(1),  and  -  a  /N*1*  . 
n  J  n  n  on 


Finally,  compute  d..  as  an  estimate  of 


Ojj  -  v.rta^}  -  -  r^2)co»{o)l,S^2)} 


according  to 


e  -  s  -  i 

°12  *11  rn  *12 


,<».«)  +  ea>e«).22 . 
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vhara  a^,  a^  *»d  a^  ••  bafora,  and 


i22  -  (n-1)'1  £(N<1)-N<1))(N<2)-N<2) 


)  . 


Than  aatloata  o  according  to 


d*  - 


e. 


a. 


2d 


12 


(N^  )2  (N*"  )2  N^N2^ 


n  n 


In  an  analogoua  manner,  we  aatloata  tha  varlanca  constant  o  (x) 
appearing  In  Equation  (10.2.12)  which  leada  to  a  confidence  Interval  for 
P{R^£x}-P{R^2^x}  .  Baaed  on  n  super-cyclaa ,  for  1*1,2  compute  d^(x) 
as  an  estimate  of 


°ll(x)  “  var{Yki>J  -  IPtR^Wcov^15,!^0} 

+  (P{R(i)5x})2var{n^i)} 


according  to 


where 


8ii(x) 


•u  <*) 


•u  ' W 


■’(%)■ 


Y(i)\2 

U<x>  +(4iy  ‘22 

1  nn  ' 


(n-1)'1  £(Y(i)-Y 

j-1  J 


.vUK2 

n  ' 
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-  n"1  S*!1*  and  n(1)  -  n_l  £]n;: 

ft  3  °  ft  3 


Finally,  compute  d.,(x)  as  an  astimata  of 


Ou(x)  -  covtY*1* ,Y^2) }  -  PiR^WcovfY^.n^} 
-  PCR^SxJcevtY^.n*2*} 

+  P{R(1)Sx}p{R(2)Sx}cov{n^1),n^2)} 


according  to 


$12 00  ■  sii^*)  “ 


/-(1)\  /— (2)\ 

■)  (x)  f  n  |a^  (a)  -1-2—1 

.W  \-(l)ri2W  [-(2)1 
\nn  /  \nn  / 

'?(1)_y(2)\ 

-(l)-(2)  ]  a22(x)  • 

®n  nn  / 


Sji^ (x) 


where  s^(x),  s2^00»  “d  *12^*^  8r8  48  befor®»  and 


a22(x)  -  (n-1)”1  C»J1> -*Jl>  >  C»^2)-«i2> ) 


Than  estimate  c  (x)  according  to 


*2,  *  m  ail(x)  +  d22(x)  2e22(x) 


5 


APPENDIX  4.  GENERATION  OF  MARKOV  CHAIN  IN  DL/I  COMPONENT  MODEL 

For  isi£7,  denote  by  X^  the  rate  perameter  of  the  exponentielly 
distributed  service  time  for  jobs  of  class  i.  Complete  knowledge  of  the 
^-matrix  of  Infinitesimal  transition  parameters  for  the  continuous  time 
Markov  chain  £»{X(t):tiO}  defined  by  Equations  (5.2.1)  and  (5.2.2)  is 
contained  in  the  routing  matrix  P  and  the  X^.  We  give  an  algorithm  for 
direct  generation  of  this  process. 

1.  Fix  an  initial  state  in  A^  and  set  t«*0. 

2.  Determine  the  status  of  the  a  and  8  center  servers,  i.e.,  whether 
or  not  they  are  busy,  and  if  so,  what  classes  of  jobs  are  in 
service.  The  8  center  server  is  busy  if  Q(t)>0;  the  a  center 
server  is  busy  if  S(t)>0. 

Assume  both  are  busy,  with  a  class  i>0  job  in  service  at  the  a  center. 

The  cases  in  which  only  one  server  is  busy  are  handled  similarly. 

3.  Generate  a  holding  time  T,  exponentially  distributed  with  rate 
parameter  X^+X^  and  advance  t  to  t'-t+T. 

4.  Determine  according  to  a  Bernoulli  trial  the  service  which 
completes  first: 

*i 

P{a  center  service  completes  first}  -  t  . 

A1  Ai 
end 

X2 

P{8  center  service  completes  first)  •  r  — r  . 

Al+ Ai 

If  a  center  service  completes  first,  go  to  6. 
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S.  Sat 

Q(t’)  -  Q(t)-1  , 

c5(f)  -  c5(t)+i  , 

S(t’)  -  S(t)  , 
and  for  2Sj*7  and  j*i,  ant 

cj(tf)  -  ^(t)  . 

If  S(t')5*7,  go  to  7.  Otharwlsa,  sat 
C?(t')  »  C?(t)+1  • 

c5(t’)  -  c5(t’)-i  , 

S(t')  -  5  , 
and  go  to  7  . 

6.  For  2£k£7,  sat 

ck(f)  -  ck(t)  . 

Generate  ie{l,2,. . .  ,7}  according  to  the  probabilities  p^,  where 
S(t)»j,  and  set 

^(t*)  -  c1(t,)+i,*i>l 

and 

Q(t')  -  Q(t)+1,  i-1  . 

Set 

S(t')  -  kQ 

and 

C.  (f)  -  C  (t')-l  , 

*0  *0 

where  kg«min{k:Ck(t ’ )>0}  . 


7.  Determine  N(t’)  as  follows.  If  the  job  completing  service  is 

the  marked  job,  set  N(t’)  equal  to  its  position  in  the  job  stack 


after  completion  of  the  service.  If  the  job  completing  service 
is  not  the  marked  job,  set  N(t')"N(t)  if  the  job  completing 
service  goes  back  in  the  job  stack  either  above  [respectively 
below]  the  marked  job,  when  prior  to  completion  the  job  was  above 
[respectively  below]  the  marked  job.  Set  N(t')*N(t)-l  if  job 
completing  service  goes  from  above  the  marked  job  to  below. 
Otherwise,  set  N(t')«N(t)+l. 

Return  to  2.  and  iterate  with  t'  playing  the  role  of  t. 
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